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ON  THE  ESTIMATION  OF  A  MOVING  SHIP’S  VELOCITY  AND  HULL 
GEOMETRY  INFORMATION  FROM  ITS  WAVE  SPECTRA 


by 

Zhijian  Wu 


Chairperson:  Guy  A.  Meadows 

The  wake  generated  by  a  moving  ship  may  extend  for  many  tens  of  kilometers  in  the 
open  ocean,  and  can  be  remotely  sensed.  Through  indirect  methods,  the  detection 
of  a  ship  and  its  related  characteristics,  is  generally  obtained  by  measuring  the  ship 
generated  waves  or  their  spectra.  From  the  viewpoint  of  remote  sensing,  interesting 
problems  exist  related  to  the  detection  of  a  ship’s  presence  and  the  acquisition  of 
dynamic  and  static  information  about  it.  This  problem  can  be  divided  into  two 
basic  aspects.  First,  how  to  obtain  a  moving  ship’s  wave  spectra  from  remotely 
sensed  images,  and  second,  how  to  extract  the  desired  ship  information  from  the 
imaged  wave  spectra.  This  thesis  concentrates  on  the  latter  aspect,  in  particular, 
how  to  estimate  a  moving  ship’s  direction,  speed,  length  and  hull  shape  from  its  wave 
spectra. 

The  extraction  of  ship  information  is  based  on  the  relations  of  the  ship’s  wave 
spectra,  wave  amplitude  function  and  hull  geometry.  In  this  thesis,  an  analytic  rep- 
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resen  tation  of  wave  elevation  is  introduced  with  the  use  of  the  Hilbert  transform, 
and  the  derivation  is  given  for  the  calculation  of  the  wave  amplitude  function  from 
the  Fourier  spectrum  of  one  and  two  dimensional  complex-valued  wave  elevations. 
Methods  and  formulas  are  given  for  estimating  a  ship’s  speed  and  direction  from  the 
spectrum  of  a  two- dimensioned  wave  patch,  a  single  wave  cut  or  two  wave  cuts.  A 
theoretical  model  of  the  wave  amplitude  function  is  developed,  and  three  methods 
are  designed  for  the  estimation  of  a  ship’s  length  from  the  wave  amplitude  func¬ 
tion.  Under  the  assumption  of  thin-ship  theory,  an  inversion  technique  to  predict 
the  geometry  of  a  ship’s  hull  from  the  wave  amplitude  function  or  its  magnitude  is  de¬ 
veloped  through  the  application  of  a  spectral  method  and  the  constrained  maximum 
likelihood  method.  Examples  comparing  theoretically  calculated  data  and  tow  tank 
experimental  data  are  given  to  demonstrate  the  methods  developed  and  estimate 
performance. 
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CHAPTER  I 


INTRODUCTION 


Although  a  ship’s  length  is  bounded  within  a  range  of  values,  the  wake  it  gener¬ 
ates  in  the  open  ocean  may  extend  for  many  tens  of  kilometers.  Through  indirect 
methods,  the  detection  of  a  ship  and  its  related  characteristics,  is  generally  obtained 
by  measuring  the  ship  generated  waves  or  their  spectra.  From  the  viewpoint  of  re¬ 
mote  sensing,  interesting  problems  exist  related  to  the  detection  of  a  ship’s  presence 
and  the  acquisition  of  dynamic  and  static  information  about  it.  For  example,  a 
ship’s  velocity,  size  and  hull  shape  are  desired  characteristics.  Figure  1.1  illustrates 
a  scheme  for  obtaining  this  information  about  a  ship  from  remotely  sensed  images. 
This  problem  can  be  divided  into  two  basic  aspects,  one  is  how  to  obtain  a  moving 
ship’s  wave  spectra  from  remotely  sensed  images,  the  second  is  how  to  extract  the 
desired  ship  information  from  the  wave  spectra.  This  work  concentrates  on  the  latter 
aspect,  in  particular,  how  to  estimate  a  moving  ship’s  direction,  speed,  length  and 
hull  shape  from  its  wave  spectra. 

Ship  generated  surface  gravity  wave  patterns  can  be  remotely  detected  using 
several  techniques,  including  visible  photography,  infrared  sensing  and  microwave 
radars.  Several  radar  remote  sensing  techniques  can  be  used  to  estimate  ocean  surface 
directional  wave  properties,  for  example,  the  Ocean  Wave  Spectrometer  (ROWS),  the 
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Figure  1.1:  A  scheme  to  extract  ship  information  from  remote  sensing  images. 

Three-Frequency  Airborne  Radar  (TRIFAR)  and  Synthetic  Aperture  Radar  (SAR) 
[1].  Vesecky  et  al.  have  studied  remote  sensing  of  ocean  waveheight  spectrum  using 
SAR  [2]- [4],  and  Monaldo  et  al.  have  studied  the  transformation  of  surface  wave 
slope-  and  height-variance  spectra  from  radar  images  [5]  in  recent  years.  These 
techniques  may  find  their  applications  in  the  detection  of  ship  characteristics  from 
remotely  sensed  images.  Ship  wave  spectra  have  distinct  features  which  are  different 
from  those  of  ambient  ocean  waves.  Tuck  et  al.  have  studied  the  Fourier  spectra  of 
real  ship  wave  elevations  and  indicated  the  possibility  of  utilizing  this  information  to 
estimate  ship  speed  [6].  The  estimation  of  ship  hull  geometry  information  from  ship 
waves  or  their  spectra  can  be  considered  as  an  inverse  Kelvin  wake  problem,  and  has 
been  explored  by  Kuhn,  Newman  et  al.  recently  [7]-[9]. 

Since  the  inverse  Kelvin  wake  problem  is  relatively  new,  little  directly  related 
published  work  exists.  The  study  of  this  problem  will  naturally  rely,  to  a  great 
extent,  on  the  existing  theories  and  results  on  the  forward  Kelvin  wake  problem,  i.e., 
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predicting  the  wave  field  due  to  the  geometry  of  the  ship.  The  forward  problem  has 
been  studied  for  several  decades  with  the  desire  to  predict  the  ship  wake  and  reduce 
the  ship  wave  resistance  to  a  ship  [15}-[28]. 

Mathematically,  the  inverse  Kelvin  wake  problem  has  some  characteristics  similar 
to  inverse  problems  in  other  areas.  Thus,  the  methods  developed  in  those  areas  will 
be  helpful  in  finding  solutions  of  this  inverse  Kelvin  wake  problem  [46]- [54]. 

In  this  work,  the  study  of  the  estimation  of  a  moving  ship’s  speed,  direction  and 
hull  geometry  characteristics  from  its  wave  spectra  is  based  on  the  recognition  of 
the  relations  among  the  wave  spectra,  the  wave  amplitude  function  and  the  ship 
hull  shape.  On  the  spectrum  diagram  the  wave  number  distribution  contains  the 
ship’s  speed  information,  the  position  of  spectrum  loci  contains  the  ship’s  direction 
information,  and  the  magnitude  of  the  wave  spectra  contains  the  ship’s  hull  geometry 
information. 

In  the  following,  Chapter  2  and  Chapter  3  serve  as  the  theoretical  foundation  of 
the  study.  Based  on  basic  ship  wave  theory  and  Fourier  theory,  Chapter  2  discusses 
ship  generated  free  waves  and  their  Fourier  spectra.  An  analytic  representation  of 
wave  elevation  is  introduced  with  the  use  of  the  Hilbert  transform.  The  concept  of 
complex  wave  elevation  with  this  analytic  representation  is  helpful  in  simplifying  the 
derivation  of  of  the  wave  amplitude  function  from  the  wave  spectra. 

In  Chapter  3,  the  derivation  is  given  for  the  estimation  of  the  wave  amplitude 
function  from  the  spectra  of  one  and  two  dimensioned  complex- valued  wave  eleva¬ 
tions.  The  spectrum  loci  on  a  spectrum  diagram,  a  distinct  characteristic  of  ship 
wave  spectra,  are  important  for  the  estimation  of  the  ship  speed  and  direction,  and 
thus  are  mathematically  described.  The  formulas  to  calculate  the  wave  amplitude 
function  are  given  in  detail.  The  effects  of  sampling  intervals  on  the  resolution  of 
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wave  angles  and  the  wave  amplitude  function,  and  the  relation  between  the  sampling 
intervals  for  a  real  ship  and  its  model  must  be  understood  and  considered  in  tow  tank 
experimental  studies  or  other  practical  applications.  Thus,  one  section  is  included 
for  these  contents. 

With  the  fundamentals  and  formulas  given  in  Chapter  2  and  Chapter  3,  the 
following  chapters  develop  the  theory  and  methodology  for  the  estimation  of  a  ship’s 
velocity,  length  and  hull  surface  shape.  Chapter  4  first  shows  the  discovery  of  the 
presence  of  a  moving  ship  in  an  ambient  random  wave  field  through  a  simple  example. 
The  methods  and  formulas  for  calculating  a  ship’s  speed  and  direction  from  the 
spectrum  of  a  two-dimensional  wave  patch,  a  single  wave  cut  or  two  wave  cuts  are 
then  discussed  in  detail.  The  examples  with  theoretically  calculated  data  and  tow 
tank  experimental  data  are  given  to  demonstrate  the  methods  and  the  estimation 
performance. 

In  Chapter  5,  a  theoretical  model  of  the  wave  amplitude  function  is  developed 
for  the  estimation  of  a  ship’s  length  from  the  wave  amplitude  function.  This  model 
explicitly  reveals  the  periodic  character  inherent  in  the  real  and  imaginary  part  and 
even  in  the  magnitude  of  the  wave  amplitude  function.  It  also  shows  the  relation  be¬ 
tween  a  ship’s  length  and  the  periodicity,  and  the  effects  of  the  bow  and  stem’s  shape 
on  the  periodic  character.  With  this  understanding,  three  methods,  the  spectrum 
method,  zero-crossing  method  and  frequency  demodulation  method,  are  designed 
to  estimate  ship  length.  Examples  are  given  for  each  method  and  the  results  are 
compared. 

Chapter  6  develops  a  technique  to  extract  a  ship’s  hull  geometry  shape  from 
the  wave  amplitude  function  or  from  its  magnitude  under  the  assumption  of  the 
thin-ship  theory.  The  spectral  method  is  used  in  converting  the  continuous  inverse 
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problem  to  a  discrete  problem,  and  the  selection  of  basis  functions  axe  discussed. 
The  ill-condition  of  matrices  in  the  resulting  equations,  noise  in  input  parameters 
and  the  wave  amplitude  function  have  severe  effects  on  the  solution  as  analyzed  in  the 
chapter.  To  reduce  the  effects  to  a  minimum,  Bayes  estimation  theory  is  applied  to 
the  inverse  problem,  and  the  constraints  are  considered  in  both  linear  and  nonlinear 
cases.  The  maximum  likelihood  method  with  constraints  is  found  to  be  especially 
useful  in  the  examples  of  mathematically  well-defined  hulls  and  that  of  a  sea-going 
tug,  the  USS  Quapaw. 

The  final  chapter,  Chapter  7,  summarizes  the  research  conducted  in  this  study 
and  gives  some  recommendations  for  further  efforts. 


CHAPTER  II 


SHIP  WAVE  SPECTRA 


The  extraction  of  ship  information  is  based  on  the  relations  of  the  ship’s  wave 
spectra,  wave  amplitude  function  and  hull  geometry.  The  wave  spectra  have  rela¬ 
tions  with  wave  elevation  and  slopes  through  Fourier  transform.  Ship  wave  numbers 
contain  the  desired  ship  speed  information,  the  position  of  spectrum  loci  in  spatial 
frequency  space  or  wavenumber  space  contains  the  ship  direction  information,  and 
the  magnitude  of  wave  spectra  contains  the  hull  geometry  information. 

In  this  chapter,  an  analytic  representation  of  wave  elevation  is  introduced  to 
simplify  the  mathematical  manipulation  in  wave  spectra  later,  and  then  the  relation 
between  elevation  and  slope  spectra  is  discussed  for  stationary  ship  wave  motions. 

2.1  Ship  Wave  Elevation  and  its  Analytic  Representation 

Propagating  waves  and  the  signals  they  carry  can  be  modeled  as  functions  of 
space  and  time,  and  they  can  be  analyzed  by  using  multidimensional  Fourier  trans¬ 
form  methods.  For  general  cases,  if  s(x,  t)  represents  a  signal  that  is  a  function  of 
spatial  position  x  =  (x,y,  z)  and  time  t,  and  S(k,u>)  represents  its  four-dimensioned 
wavenumber  frequency  spectrum,  then  s(x,<)  and  5(k,w)  can  be  expressed  in  terms 
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Figure  2.1:  Reference  coordinate  system.  The  mean  water  surface  is  at  z  —  0. 

of  each  other  in  the  following  equations  [10]  : 

S(k,a;)=  f°°  f°°  s(x,t)e-j(u,-kx)dxdt  (2.1) 

s(x,t)  =  r  r  (u,-lc'x) dkdu)  (2.2) 

47T  7—oo  7—oo 

where  j  =  -y/— I,  and  k  •  x  represents  the  inner  product  of  the  wavenumber  k  and 
the  position  vector  x.  The  space-time  signal  s(x,  t)  can  be  considered  as  the  su¬ 
perposition  of  numerous  elemental  propagating  waves  exp  {j(u>  —  k  •  x)  weighted  by 

S(k,w). 

In  the  problems  below,  the  wave  field  is  generated  by  a  moving  ship  in  deep  water 
and  it  can  be  described  as  a  three-dimensional  problem,  i.e,  x  and  y  in  space  and  t 
in  time.  Additionally,  the  assumption  of  linearized  free-surface  boundary  condition 
is  made.  Now,  consider  a  reference  system  moving  with  the  ship  in  the  positive 
x-direction  with  speed  U ,  as  shown  in  Figure  2.1,  then  the  wave  elevation  r)(x,y,t) 
can  be  expressed  in  the  following  form  [14] 

r)(x,y,t)  =  Re{ J^dw  deA(^e)e-^XCM0+v^e)+^-KUcMe)^  .  (2.3) 
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Here,  0  is  the  wave  angle,  A{u,0)  is  called  the  wave  amplitude  function,  and  K{uj) 
is  the  wave  number  corresponding  to  a  given  frequency  u  in  accordance  with  the 
dispersion  relation  for  infinite  depth 

K  =  7  <2'4> 

where  g  denotes  the  acceleration  of  gravity.  For  real  problems,  the  wave  elevation  is 
always  real;  thus,  the  operation,  to  take  the  real  part  is  used  in  (2.3). 

It  will  be  more  convenient,  however,  if  the  real  operation  in  (2.3)  can  be  left  out  in 
the  complicated  mathematical  manipulation.  For  this  purpose,  the  complex- valued 
wave  elevation  is  introduced  here  with  the  use  of  the  Hilbert  transform.  The  Hilbert 
transform  of  a  real-valued  signal  x(^)  is  another  real-valued  signal,  which  is  defined 
by  the  convolution  integral  of  x(i)  and  [11].  That  is,  if  the  Hilbert  transform  of 
x(t)  is  denoted  by  x(t),  then 

#*)-*{  *<«»  =  £^7)*-  <2'5> 

An  analytic  signal  x(t)  associated  with  x(t)  is  defined  by 

x(0  =  x(0  +  j  X(0  (2-6) 

and  it  can  be  expressed  with  a  magnitude  function  a(t )  and  a  phase  function 
where  a(t)  describes  the  envelope  of  the  original  function  x(*)>  and  y?(<)  describes 
the  instantaneous  phase  of  x(0-  Thus,  x(f)  can  be  written  in  the  form 

m  =  «w  (2.7) 


a(0  =  [x(<)  +  x(0J* 

<f(t)  =  tan-1[  ]  . 


where 
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A  useful  property  of  an  analytic  signal  x(t)  is  the  Fourier  transform  relations 
among  x(0»  x(*)  and  x(0-  Let  X(f),  X{f)  and  X(f)  be  the  Fourier  transforms  of 
x(<),  x(f)  and  x(<).  respectively.  Then,  it  can  be  proved  that 


X(f)  =  -j  sgn(f)X(f) 

► 

(2.8) 

2  X(f) 

for  /  >  0 

X(f)  =  [1  +  sgn(f)]X(f)=<  X{f) 

for  /  =  0 

(2.9) 

0 

. 

for  /  <  0 

where  sgn(f)  is  a  sign  function  which  is  defined  as  sgn(f)  =  1  if  /  >  0,  sgn(f)  =  0 
if  /  =  0  and  sgn(f)  =  — 1  if  /  <  0.  From  these  relations,  X(f)  and  X(f)  can  be 
obtained  once  X(f)  ic  available. 

Now,  consider  two  analytic  signals  associated  with  the  real  signals 

/i(t)  =  i2e{Ae-jnT}  (2.10) 

Mr)  =  Re{  Ae^  )  (2.11) 

where  A  is  complex  and  independent  of  t,  and  fl  is  real.  With  the  help  of  the 
following  Hilbert  transform  relations  of  sine  and  cosine  with  constants  c\  and  c2: 

H{  cos(cjT  +  c2) }  =  sin(cir  +  c2) 

H{  sin(ciT  +  c2)  }  =  —  cos(cir  +  c2) 

the  analytic  signals  of  fi{r)  and  /2(r)  axe  given  by 

A*ejnT  for  n  >  0 

/i(T)  =  "  (2.12) 

Ae~jnT  for  <  0 


/2(r)  = 


Ae^Ur  for  ft  >  0 
for  ft  <  0 


(2.13) 
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where  A *  denotes  the  conjugate  of  A. 

By  applying  the  above  concept  for  an  analytic  signal  to  the  wave  elevation  given 
in  (2.3),  the  complex- valued  wave  elevation,  i.e.,  the  analytic  representation  of  the 
wave  elevation  can  be  derived.  In  terms  of  formula  (2.12),  the  analytic  representation 
of  rj(x,y,t)  is  now  given  by 

y{x,y,t)  =  T)(x,y,t)  +  iii(x,y,t) 

=  J°°<L>  J*  de  Am(u,e)eiK{xcoa6+yaine)~i{u,~KUco,l9)t 

+J°°du  d6  A{u,0)e-iK(xc~9+v™6)+i('d'-KUco'6)t  (2.14) 

where  the  Hilbert  transform  is  taken  with  respect  to  x.  The  above  analytic  repre¬ 
sentation  of  wave  elevation  and  the  property  of  the  Hilbert  Transform  in  (2.9)  will 
be  helpful  in  the  following  mathematical  derivation. 

2.2  Spectra  of  Wave  Elevation  and  Slope 

Equation  (2.3)  describes  a  non-steady  wave  motion,  that  is,  the  wave  elevation 
changes  not  only  with  the  spatial  position  but  also  with  time.  If  the  motion  is  steady 
state  in  the  reference  system  moving  with  the  ship,  however,  expression  (2.3)  must 
be  independent  of  time.  Thus, 

KU  cos0-u  =  0  (2.15) 

and  the  phase  velocity  of  each  admissible  wave  component  can  be  obtained  from 
(2.15), 

V,  =  ^-=U cos#  . 

K 


(2.16) 
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From  (2.4)  and  (2.15),  the  wave  number  K(9 )  can  be  related  to  the  ship  speed  and 
wave  angle  9  in  the  form 

"  whro  ■  <2'17> 

The  restriction  (2.15)  can  be  used  to  eliminate  the  variable  t.  By  retaining  the  wave 
angle  0  and  by  noting  that  (2.15)  and  (2.16)  require  that  cos  6  >  0,  the  free-wave 
distribution  of  a  given  ship  for  the  deep-water  case  now  can  be  expressed  from  (2.3) 
in  its  real  form 

ri(x,y)  =  Re{f*  A(0)c-iVCa<0>~l’K‘<0Md0}  (2.18) 

where 

Kx(0)  =  K{0)cos9 
Ky{9)  =  K{0)  sinfl. 

Note  that  0  must  now  range  from  — ^  to  *  because  of  the  requirement  cos  0  >  0. 
The  analytic  representation  of  wave  elevation  now  becomes 

rj(x,y)  =  J*wA*(O)ei[K‘Wx+K’Wv]d0  (2.19) 

where  the  Hilbert  transform  is  taken  with  respect  to  x.  Additionally,  the  wave 
elevation  spectrum  is  given  for  the  real  wave  elevation  by 

H(u,  v)  =  JF{  :,(*,  y) }  =  r  r  n{x,  y)e-^+^dxdy  (2.20) 

J — oo  J— 00 

and  for  the  complex  wave  elevation  by 

H( u,v)  =  *■{$(!,»)  }  =  r  r  Hx,y)e-M**-'dxdy  (2.21) 

J—oo  J—oo 

where  u  and  v  are  the  spatial  frequencies  associated  with  x  and  y,  respectively. 
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The  relation  between  the  wave  elevation  and  its  spectrum  has  been  now  estab¬ 
lished  through  Fourier  relations.  Theoretically,  thus,  the  ship  information  can  be 
recovered  from  either  the  wave  elevation  itself  or  its  spectrum.  In  some  real  situ¬ 
ations  such  as  in  remote  sensing,  however,  the  available  information  on  ship  waves 
may  not  be  the  wave  elevation  or  its  spectra  but  wave  slopes  or  the  slope  spectra. 
Therefore,  a  brief  investigation  of  the  relations  among  them  is  made  below. 

The  surface  wave  slope  vector  s  is  defined  by  the  gradient  of  the  wave  elevation, 
i.e., 

s(x,y)  =  Vj;  -  Tjx(x,y)x  +  fjy(x,y)y  (2.22) 

where  V  denotes  the  gradient  operation,  x  and  y  denote  the  unit  vector  in  x—  and 
y—  directions,  and  fjx  and  fjy  are  the  partial  derivatives  of  wave  elevation  with  respect 
to  x  and  y,  respectively.  In  terms  of  the  properties  of  Fourier  transform,  the  spectra 
of  the  slope  components  are  given  by 

Hx{u,v)  =  Jc{  fjx(x,  y)  }  =  j2jrui/(u,  v)  (2.23) 

Hy(u,v )  =  F{fjy{x,y)}  =  j2irvH(u,v)  (2.24) 

and  the  vector  slope  spectrum  is  given  by 

S(u,u)  =  .F{s(x,y)  }  =  j27r(ux  +  vy)H(u,v)  .  (2.25) 

A  slope  component  in  any  direction  n,  denoted  by  sn(x,  y),  can  be  obtained  from  the 
directional  derivatives  of  wave  elevation,  that  is, 

5„(x,  y)  =  n  •  S(s,  y)  .  (2.26) 


Thus,  its  spectrum  is 

Sn(u,  v)  =  n  •  !F{  s(x,  y)  }  =  j27r  [n  •  (ux  +  vy)]H(u,v)  . 


(2.27) 
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In  terms  of  the  above  relations,  slope  spectra  can  be  obtained  once  elevation  spectra 
are  available.  Similarly,  elevation  spectra  is  also  obtainable  from  slope  spectra,  except 
at  zero  frequency  since  the  denominator  will  be  zero  at  zero  frequency  .  The  zero 
frequency  component  of  the  elevation  spectra  represents  the  average  value  of  wave 
elevation,  and  it  is  not  very  important  in  many  real  situations.  Because  of  the 
relation  between  slope  spectra  and  elevation  spectra,  the  discussion  in  later  chapters 
will  focus  only  on  the  elevation  spectra. 


CHAPTER  III 


ESTIMATION  OF  THE  WAVE  AMPLITUDE 
FUNCTION  FROM  WAVE  SPECTRA 


In  this  chapter,  ship  wave  spectrum  patterns  are  discussed  and  the  relationship 
between  the  wave  amplitude  functions  and  wave  spectra  is  derived.  As  stated  in 
the  last  chapter,  the  ship  generated  stationary  wave  elevation  and  its  spectrum  are 
entirely  dependent  on  wave  numbers  and  the  wave  amplitude  function.  Generally 
speaking,  the  wave  elevation  or  wave  spectra  can  be  directly  measured  or  remotely 
sensed,  but  the  wave  amplitude  function  can  not.  As  will  be  seen  in  the  following 
chapters,  however,  dynamic  and  static  information  about  a  moving  ship,  such  as 
speed  and  hull  geometry,  is  strongly  reflected  in  wave  numbers  and  the  wave  ampli¬ 
tude  function.  Thus,  recovering  the  wave  amplitude  function  from  wave  spectra  is 
an  important  procedure  for  obtaining  this  information. 

The  wave  pattern  analysis,  especially  the  estimation  of  the  wave  amplitude  func¬ 
tion,  also  plays  an  important  role  in  the  ship  wave  resistance  analysis  because  of  its 
direct  relation  with  wave  resistances.  For  the  purpose  of  analyzing  the  ship  wave 
resistances,  different  derivation  methods  to  calculate  the  wave  amplitude  function 
were  introduced  in  the  last  twenty  years  [15]-[22].  Ship  wave  patterns  not  only  con¬ 
tain  the  information  that  can  be  extracted  to  estimate  ship  wave  resistances,  but 
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also  contain  other  useful  ship  information,  thus,  ship  wave  pattern  analysis  is  also 
tin  important  means  for  remotely  sensing  ship  information  [6]. 

In  this  chapter,  the  explicit  and  succinct  expressions  for  the  wave  amplitude 
function  are  established  from  one-  or  two-dimensional  wave  spectra  based  on  the 
analytic  representation  of  wave  elevation.  The  distinct  characteristics  of  ship  wave 
spectra  can  be  observed  in  these  derivations,  and  they  become  the  basis  for  the 
estimation  of  ship  speed  and  direction  in  the  next  chapter.  The  effect  of  sampling 
intervals  on  the  wave  angle  and  wave  amplitude  function  is  also  discussed  for  practical 
use  of  this  formulas.  In  this  chapter,  the  first  section  discusses  two  dimensional  wave 
fields;  the  second  section  discusses  one  dimensional  wave  cuts;  and  the  final  section 
discusses  the  effect  of  sample  intervals. 

Before  discussing  the  estimation  of  the  wave  amplitude  function,  it  is  helpful  to 
review  briefly  some  formulas  about  the  6  function,  which  can  be  found  in  reference 

[41], 


*(*) 

=  6(-x) 

(3.1) 

S(x,y) 

=  *(*)%) 

(3.2) 

f(x)S{x  -  o)  =  f(a)6(x  -  a)  (3.3) 

J  S(x  -  y)6(y  -  a)dy  =  S(x  -  a)  (3.4) 

where  Xj  is  the  root  of  q(x)  =  0,  and  its  derivative  q'(xj)  ^  0  . 


(3.5) 
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3.1  Estimation  from  2-D  Wave  Fields 

As  seen  in  Chapter  2,  the  ship  generated  surface  wave  can  be  described  as  a 
two-dimensional  (2-D)  wave  field.  The  two-dimensional  wave  or  spectrum  data  can 
be  obtained  either  in  tow  tank  experiments  [23],  [30],  or  possibly  by  remote  sensing 
methods  [5],  [31].  The  most  distinct  signature  contained  in  ship  wave  spectra,  dif¬ 
ferent  from  general  ocean  wave  patterns,  is  the  locus  in  the  spatial  frequency  plane. 
The  first  subsection  below  discusses  the  properties  of  the  wave  spectra.  The  second 
subsection  discusses  the  relation  between  the  wave  amplitude  function  and  the  wave 
spectra,  the  discrete  forms  for  calculating  the  wave  amplitude  function  and  the  effect 
of  truncated  errors. 

3.1.1  Loci  of  Ship  Wave  Spectra 

In  this  subsection,  discussion  starts  from  the  complex-valued  wave  elevation. 
From  Chapter  2,  a  stationary  wave  elevation  is  given  by 

rj{x,y)  =  »?(x,y)  +  j  fi{x,y)  =  / 2  Am{e)e^Kz(B)xJrKx(6)v]d9  (3.6) 

where  the  wave  angle  6  ranges  from  —  |  to  ^  and  the  wave  amplitude  function  A(6) 
is  complex. 

To  obtain  the  wave  elevation  spectrum,  the  Fourier  transform  of  the  wave  eleva¬ 
tion  is  taken  and  it  follows  that 

H(u,v)  =  r  r  r}(X,y)e-™^dxdy 

J—oo  J— oo 

=  j\  Am(0){ J°°  f  e~i2*Uu-]£iPMv-!%r1)y]dxdy}d0 
=  A*(6)6(u  -  -  ^£p-)d0  . 


(3.7) 
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Now,  let 


gx{9)  =  u- 


g3(0)  =  v  - 


KM 


KM 


and  let  9\  and  92  be  the  roots  of  gi(0)  and  <72(0),  respectively.  Then,  in  terms  of  the 
8  formulas,  (3.7)  becomes 

fr.  . *(»-*)*( 

-  wJmm6(e'-h)  (3-10> 

where  <7X  and  g2  represent  the  derivatives  of  gx  and  g2,  respectively,  with  respect  to 


(3.10) 


9.  In  terms  of  (3.10),  thus,  the  wave  spectrum  H(u,v)  is  combined  with  a  number 
of  impulses  with  intensity  ,  1,*  at  position  9x  =  02.  Here,  9X  or  02  is  solved 

from  <7i(0)  =  0  and  g2(9)  =  0,  that  is, 


«  =  Ml 

2tt 

v  =  Ml 
2?r 


(3.11) 


(3.12) 


This  set  of  parametric  equations  represent  a  locus  on  the  (u,  v)-plane,  or  the  wave 
spatial  frequency  plane  and  describe  the  distribution  of  ship  wave  components  on 
this  plane.  That  is,  only  the  spectrum  values  on  this  locus  are  non-zero.  Besides, 
the  spatial  frequencies  u  and  v  in  the  domain  of  the  Fourier  transform  are  consistent 
very  well  with  the  ship  wave  numbers,  Kx  and  Ky,  in  the  z—  and  y—  directions. 
Thus,  the  locus  on  the  spatial  frequency  plane  contains  ship  information. 

The  wave  angle  9  and  the  wave  number  K(0)  can  be  expressed  as  functions  of  u 
and  v.  By  solving  (3.11)  and  (3.12),  it  follows  that 


9  =  tan  *(-) 
u 


(3.13) 


X  I — . — . — . — i—  -A... _ , A _ . _ i _ _ _ _ _ _ _ _ _ , _ 2^ _ _ _ _ _ I 

0.0  0.1  0.2 


u  (1/m) 

Figure  3.1:  Loci  of  wave  spectra  in  the  spatial  frequency  (u,  v)  plane  for  different 
ship  speeds  U. 

K(9)  =  -^  =  2ir(u2  +  u2)*  (3.14) 

cos*  a 

where 

=  ■  (3-15) 

The  locus  of  ship  wave  elements  can  be  also  expressed  by  one  equation.  For  this, 
canceling  0  in  (3.13)  and  (3.14),  it  follows  that 

+  =  °  .  (3.16) 

The  loci  on  the  (u,  v)  plane  for  different  moving  speeds  are  plotted  in  Figure  3.1. 


As  shown  in  the  figure,  the  wave  components  are  located  in  the  first  and  fourth 
quadrants  of  the  (u,  v)  wave  spatial  frequency  plane  since  the  wave  angles  9  are  in 
(— f ,  f]  and  u  >  0. 
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3.1.2  Calculation  of  the  Wave  Amplitude  Function 

In  this  subsection,  an  explicit  and  succinct  expression  for  the  wave  amplitude 
function  is  established  from  two-dimensional  continuous  or  discrete  wave  spectra. 

Since  the  wave  spectrum  is  combined  with  the  8  functions  as  shown  in  (3.10), 
it  is  possible  to  obtain  the  expression  of  A(9)  by  integrating  both  sides  of  (3.10), 
resulting  in 

A'(0,)  =  l*(»i)llft(«i)l  .  (3.17) 

However,  this  result  complicates  further  manipulation,  because  the  integral  is  with 
respect  to  02  and  the  integrand  is  a  function  of  tt  and  v,  although  there  is  the  relations 
among  u,  v  and  02. 

An  alternative  method  is  to  start  the  derivation  directly  from  (3.7)  and  integrate 
with  respect  to  u  since  H(u,v)  is  a  function  of  u  and  v.  In  terms  of  the  properties 
of  the  6  function,  (3.7)  can  be  written  as 

£(»,»)  =  A-(«„)f^S(u-^)S(v-^)de  (3.18) 

where  0O  =  0o(u,  v)  satisfies  (3.11)  and  (3.12).  In  order  to  simplify  this  expression, 
let 


m  = 


Ky(9)  sin# 

2ir  °2ircos20 


then  (3.18)  becomes 


H{u,v)  =  Am(0 0)  /*  S(u  -  R  cot  9)6{v  -  R)d0 
•M 


A-(«„)  j_ 
A"(t o) 


00  1 


r:(0o) 


K(9) 

8(u  —  vh(v)) 


6{u  -  Rh(R))8(R  -  v)dR 


(3.19) 
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where  h(R)  =  cot  9  and  the  derivative  R'  is  given  by 


Tt,(Q\  -  dRW  _  K  1  +  sip2  6 
R{-9’  d0  K°  2ir  cos3  9  ‘ 


(3.20) 


Note  that  u  —  Rh(R)  =  0  and  R  —  v  =  0  are  equivalent  to  (3.11)  and  (3.12);  thus, 
the  last  step  in  (3.19)  is  valid.  Here,  Rh(R)  is  given  by 


Rh(R)  ss  Rcot9 

=  J^(K0  +  yjKl  +  16x2/P)  .  (3.21) 

To  cancel  the  6  function  on  the  right  side  of  (3.19),  integrate  both  sides  of  (3.19) 
with  respect  to  u.  The  result  is 

£>■’>**=$$  •  (322) 
Since  H(u,v)  is  combined  with  6  functions  as  shown  in  (3.19),  it  is  possible  to  write 
H(u,  u)  into  the  6  function  with  its  intensity  Hint(u,  v),  i.e.,  H(u,  v)  =  Hint{u,  v)S{u  — 
vh(v)).  Thus,  the  integral  of  the  left  hand  side  of  (3.22)  is  equal  to  Hint(vh(v),  v ),  and 
it  represents  the  intensity  of  the  spectrum  on  the  locus.  Hence,  the  wave  amplitude 
function  can  be  written  from  (3.22)  in  the  form 


y4(0(u,v))  =  R'(9{u,v))H;nt(u,v) 

=  jr„l . 

2ircos  39(u,v)  n  7 


(3.23) 


Note  that  u  is  retained,  and  9’s  subscript  “0”  is  omitted  in  (3.23)  for  simplicity,  but 
remember  that  u  and  v  must  be  on  the  locus,  that  is,  must  satisfy  (3.11)  and  (3.12). 

The  wave  amplitude  function  is  an  even  function  of  the  wave  angle  when  the 
ship  wave  is  symmetric  in  the  y-direction.  This  can  be  proved  in  terms  of  the  above 


relations.  From  (3.13),  the  wave  angle  is  a  odd  function  of  v,  i.e.,  9(u,  v)  =  —9(u,  — v ). 
When  the  ship  wave  is  symmetric  in  the  y-direction,  that  is,  rj(x,y)  =  rj(r  — y )  the 
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wave  spectrum  is  even  with  respect  to  v,  i.e.,  H(u,v)  =  H(u,  — 1>),  according  to  the 
Fourier  transform  properties.  Thus,  H(u,  v)  =  H(u ,  —v)  too.  It  is  found  from  (3.23) 
that  the  amplitude  function  A(0)  is  even  with  respect  to  0,  that  is,  A(6)  =  A(— 6), 
when  the  wave  elevation  is  even  with  respect  to  y.  For  many  ship  types,  the  hull  is 
symmetric  with  respect  to  the  ship  central  plane,  thus  the  ship  wave  is  symmetric 
with  respect  to  the  centerline,  and  the  wave  amplitude  function  is  even  with  respect 
to  wave  angles. 

There  are  different  ways  of  obtaining  the  wave  spectrum  in  real  applications.  One 
way  is  from  the  wave  elevation  or  slopes,  that  is,  the  wave  spectrum  is  calculated 
by  taking  the  Fourier  transform  of  the  observed  wave  data.  Another  possibility  is  to 
obtain  the  wave  elevation  spectra  indirectly.  For  instance,  the  wave  spectrum  can  be 
estimated  from  radar  images.  In  practical  applications,  the  data  are  discrete;  thus, 
the  discrete  formula  is  useful  for  real  situations  and  a  discussion  is  given  below. 

To  derive  the  discrete  formulas,  some  definitions  of  discrete  variables  are  given 
first.  If  we  let  ni  and  n2  be  the  discrete  forms  of  x  and  y,  and  Ax  and  Ay  the 
interval  sizes  in  the  x—  and  y— directions  in  the  spatial  domain,  and  let  fcj  and  k2 
be  the  discrete  forms  of  u  and  v  ,  and  Au,  Av  the  interval  sizes  in  the  u—  and  in¬ 
directions  in  the  spatial  frequency  domain,  then  x  =  niAx,  y  =  n2Ay,  u  =  fcjAu, 
and  v  =  k2Av. 

Under  the  assumption  that  fj(x,  y)  is  very  small  outside  the  range  of  —  hfr  <  x  <  k*- 
and  the  infinite-range  Fourier  transform  can  be  approximated  by  its 

Fourier  transform  with  finite  ranges  L\  and  L2  in  the  x—  and  y— directions,  i.e., 

H{u,v)  =f  f  ij(x1y)e~^2^ux'¥vy^dxdy 

j— OO  J — oo 

*  /T, .  (3.24) 

2  t 
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After  discretization,  (3.24)  can  be  expressed  in  the  form 


Sx-i  Z2.-1 

H(ki,k2)  «  AxAy  £  H  »H*i>  fc2)exp[-j2ff(^p-  +  — p)]  (3.25) 

^n2=-^  iVl  ^2 


ni 


*2  = 


iV, 
2  ’ 


—1,0, 1,2,...,-—  —  1 

—10  12  —  —  l* 

AjVjllA,  **.y  ^ 


where  iVx  and  N2  are  the  number  of  data  points  in  the  x—  and  y—  directions. 

The  right  hand  side  of  (3.25)  can  be  computed  using  powerful  FFT  algorithms. 
Once  the  discrete  Fourier  transform  of  wave  elevation  has  been  found,  the  algorithm 
for  recovering  the  wave  amplitude  function  can  be  established  from  (3.23)  in  the 
discrete  form 


A(6(k  !,*,))  =  R\9(kuk2))H-(kuk2)Au 
1  +  sin2  9(ki,  k2) 


=  K0 


2xcos36(ki,k2) 


H'(  kuk2)Au 


(3.26) 


where  K0  —  -fa,  A u  —  ,  and  0 ,  k\  and  k2  satisfy 


9(kuk2)  =  tan_1(^)  (3.27) 

ki 

(A«*,)’  +  (A»*,)’  =  =  (3.28) 

That  is,  the  spectrum  values  cam  only  be  taken  from  those  on  the  spectrum  locus. 
Equation  (3.26)  does  not  contain  the  summation  operation  because  there  is  only  one 
non-zero  value  of  H(u,  v)  for  each  v.  Once  the  discrete  spatial  frequency  spectrum 
H(ki,k2)  of  the  ship  wave  elevation  fj(x,y)  is  obtained,  the  ship  wave  amplitude 
function  can  be  extracted  from  it. 

In  practical  cases,  the  real- valued  ship  wave  elevation  rj(rii,n2)  or  its  spatial 
frequency  spectrum  H(k\,k2)  is  available.  According  to  the  definition  of  77(1,3/) 


i 
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and  the  properties  of  the  Hilbert  transform,  the  spectra  H(kx,  k2)  and  H(k\,  k2)  are 
related  by 

* 

2H(ki,k2)  for  ki  >  0 

H(kuk2)=<  H(ki,k2)  forjb,=0  (3-29) 

0  for  ki  <  0 

Thus,  the  wave  amplitude  function  also  can  be  expressed  by  the  spectrum  H(ki,k2), 
that  is, 

A(«(*„Jfc,))  =  K„-  (3-30) 

7T  COS3  0{ki,  k2) 

with  fcj  >  0.  Note  that  the  value  at  k2  =  0  is  not  considered  in  the  above  formula, 
since  when  ki  =  0,  9  =  and  cos  9  =  0.  This  will  result  in  the  infinity  of 

A{9{kuk2)). 

In  the  above  discussion,  it  has  been  assumed  that  the  truncation  error  caused  by 
the  finite  data  length  in  the  x—  and  y—  directions  can  be  neglected,  and  thus  the 
FFT  algorithm  is  used  to  obtain  the  wave  spectrum,  and  then  the  wave  amplitude 
function  is  recovered  from  it.  In  some  situations,  however,  the  truncation  error  is  too 
large  to  be  neglected.  In  this  case,  the  wave  amplitude  function  may  be  recovered 
by  using  an  inversion  technique.  It  is  assumed  that  the  truncated  wave  elevation  is 
represented  by 

Vr(x,y)  =  fi{x,y)gT(x,y)  (3.31) 


Gt(u,v)  =  LiL2Sa(irLiu)Sa(nL2v) 


(3.33) 
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where  Sa(()  =  sin£/£  is  called  the  sampling  function.  By  using  the  convolution 
property  and  the  result  in  (3.7),  the  spectrum  of  the  truncated  wave  elevation  in 
(3.31)  is  given  by 

Ht( «,  v)  =  LxL2  /_*  Am(9)Sa(rLi(u  -  ^l))Sa(irL2(v  -  ^p-))dd  .  (3.34) 

Thus,  if  Ht(u,  v)  is  known,  the  wave  amplitude  A{9)  can  be  founded  by  solving  this 
integral  equation. 

3.2  Estimation  from  1-D  Wave  Cuts 

In  many  real  situations,  two-dimensional  wave  elevation  data  or  spectra  may  not 
be  available,  but  one-dimensional  (1-D)  wave  cuts  or  spectra  may  be  obtained  by 
some  remote  sensing  means  or  by  field  measurements.  In  this  section,  the  relation 
between  one-dimensional  wave  spectra  and  the  wave  amplitude  function  is  discussed 
based  on  the  complex-valued  wave  elevation  cuts.  Then,  a  special  example  is  given, 
in  which  the  wave  cut  is  measured  by  a  stationary  sensor,  and  the  wave  amplitude 
function  is  recovered  from  its  FFT  spectrum. 

The  derivation  of  the  wave  amplitude  function  from  wave  spectra  includes  three 
steps.  First,  the  one-dimensional  wave  cut  is  expressed  as  a  function  of  time  accord¬ 
ing  to  a  pair  of  wave  cut  path  equations;  then  the  spectrum  is  obtained  by  taking  the 
Fourier  transform  of  the  wave  elevation;  finally,  the  wave  amplitude  function  is  ex¬ 
pressed  in  terms  of  the  spectrum.  Now,  consider  a  general  case  shown  in  Figure  3.2, 
where  a  wave  cut  is  taken  in  the  ship  generated  wave  field  by  a  sensor  moving  with 
a  uniform  speed  Up  in  a  direction  making  an  angle  a  with  the  positive  x-axis,  i.e., 
the  direction  that  the  ship  is  moving.  In  the  reference  system  moving  with  the  ship 
at  constant  speed  U,  the  wave  cut  path  can  be  described  by  a  pair  of  parametric 
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Figure  3.2:  1-D  wave  cut  in  the  ship  wave  field. 


equations 

x(t )  =  xo  -f  (Up  cos  a  —  U)  t  (3.35) 

y(t)  =  y0  +  UpSinat  (3.36) 

where  t  is  a  parameter  representing  the  measurement  time,  and  (xo,  yo)  is  the  sensor’s 
initial  position  in  the  given  coordinate  system  at  t  =  0.  If  the  sensor  is  mounted 
on  an  airplane  or  satellite,  the  sensor’s  speed  Up  will  be  much  larger  than  the  ship’s 
traveling  speed  U.  In  some  cases,  however,  the  sensor  is  considered  to  be  fixed  in  a 
position  to  measure  the  wave  cuts  when  a  ship  passes  through,  for  instance,  in  tow 
tank  experiments.  When  Up  cos  a  >  U,  x(t)  increases  with  t,  otherwise  it  decreases 
with  t. 

The  complex-valued  wave  elevation  cut  can  be  derived  by  substituting  (3.35)  and 
(3.36)  into  (2.18)  and  then  by  taking  the  Hilbert  transform.  Its  form  is  given  by 

m  =  «(*(«).  »(0)  =  { 


(3.3T) 


26 


Figure  3.3:  Curves  of  versus  0  for  different  wave  cut  angles  a. 

where  {.}**  specifically  defines  a  conditional  conjugate  operation,  that  is,  no  opera¬ 
tion  is  taken  when  3>(0)  >  0,  but  the  conjugate  operation  is  taken  when  $(#)  <  0. 
$o(0)  and  $(0)  are  defined  by 

*o(«)  =  Kx{0)xo  +  Ky{d)y0  (3.38) 

*(6)  =  ~  [Kx{0)(Up  cos  a -U)  +  Ky(0)Up  sin  a  ] .  (3.39) 

i'K 

With  the  above  wave  cut  expression,  the  wave  spectrum  is  obtained  by  taking 
the  Fourier  transform  of  i?(£)  with  respect  to  t: 

H(f)  =  )}= 

J  —  oo 

=  J\  {/!•(<>)  }••  S(J  -  |*(»)| )M 

=  {/!■(«„)  (3,10) 
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Figure  3.4:  Curves  of  versus  0  for  different  wave  cut  angles  a. 
where  0o  must  satisfy  the  equation 


/  -  W)l  =  0  (3.41) 

and  must  be  such  that  $  (^o)  ^  0.  Here,  the  derivative,  4>’(0),  is  calculated  from 
(3.39)  and  is  given  by 

$'(0)  =  —  — —  [((/_  cos  a  —  C/)sin0  +  £/p  sin  a  sec  0(sin2  0  +  1)]  .  (3.42) 

2 iru*  cos*  0 

When  Up  cos  a  >  U,  $(0)  and  $'(0)  can  be  approximated  by 

KqUp  cos (0  —  a) 


$(0)  « 
$'(0)  « 


2  ir  cos2  0 

r 

y— [  cos  a  sin  $  +  sin  a  sec  0(sin3  9  +  1)] 


(3.43) 

2x  cos2  0 1 -  - v .  '  *,J  (3'44) 

where  Kq  =  -fa .  These  forms  are  similar  to  the  ones  given  by  Tuck,  Collins  and  Wells 

[6].  The  curves  of  0  versus  and  are  plotted  in  Figure  3.3  and  Figure  3.4 
for  the  above  approximate  relations. 
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By  rewriting  (3.40),  the  wave  amplitude  function  A{0)  now  can  be  expressed  in 
the  form 

A{9)  =  |*'(0)|  {  Hm(f)  }*♦  (3.45) 

where  the  subscript  of  0  is  omitted  for  simplicity,  but  keep  in  mind  that  9  must 
satisfy  (3.41).  When  the  ship  wave  is  symmetric  to  the  ship  central  line,  the  wave 
amplitude  function  is  an  even  function  of  9,  thus  the  calculation  may  be  taken  only 
for  0  <  0  <  j  or  <  0  <  0.  Because  of  the  fact  that  H(f)  =  0  for  /  <  0,  only 
the  positive  frequency  components  needs  to  considered  in  the  above  calculation.  The 
valid  wave  angle  can  be  found  for  each  positive  /  by  solving  (3.41)  or  by  estimating 
from  curves  in  Figure  3.3  and  Figure  3.4.  In  terms  of  (3.45),  the  wave  amplitude 
function  can  be  calculated  if  the  wave  cut  spectrum  H(f)  has  been  obtained  together 
with  the  parameters  U,  Up,  a,  and  the  initial  position  (x0,yo)-  The  methods  for 
obtaining  the  ship  speed  U  and  the  wave  cut  angle  a  will  be  presented  in  the  next 
chapter. 

As  an  example,  a  special  case  now  is  considered  in  which  a  sensor  is  assumed  fixed 
in  position  and  the  wave  elevation  is  measured  as  the  wave  field  is  generated  when 
the  ship  passes,  as  is  typically  done  in  tow  tank  wave  measurement  experiments. 
Both  Up  and  a  are  zero  for  this  case.  0  <  cos  9  <  1  when  0  ranges  from  to 
thus,  now  $(0)  =  —  i*Ucx* 9  <  0*  For  this  case,  the  wave  amplitude  function  in  (3.45) 
and  the  /  constraint  condition  in  (3.41)  become 

A{9)  =  =  ff(/)^^ej*o(0)  (3.46) 

ZirU*  cos*  v 

f _ 9— 

2irU  cos  0 

Since  frequency  /  is  positive,  its  minimum  value  that  can  be  considered  in  the  esti¬ 
mation  of  the  wave  amplitude  function  from  the  spectrum  is  /mtn  —  ^7  in  terms  of 
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(3.47).  Given  a  value  of  /  larger  than  /m,B,  there  are  corresponding  values  for  wave 
angle  B  and  //(/).  Hence,  A(B)  can  be  evaluated  in  terms  of  (3.46). 

Usually,  the  spectrum  H(f)  is  obtained  from  the  FFT  of  discrete  wave  cut  data  or 
mapped  from  other  discrete  spectra.  Therefore,  /  is  discrete.  For  the  wave  elevation 
cut,  if  the  data  length  is  N  and  the  sample  interval  is  At,  then  the  discrete  positive 
frequency  is  given  by 

f  =  k 
Jk  NAt 

k  —  IcQ)  ko  4”  1, ...,  ^ 

where  ko  denotes  the  smallest  integer  larger  than  g^-.  The  corresponding  discrete 
wave  angle  is  given  by 


*  -  “*"0 


k  =  ko,ko  +  1,...,—  . 


(3.49) 


From  this  formula,  it  is  found  that  the  data  length  should  be  large  in  order  to  obtain 
a  good  resolution  for  small  wave  angles.  The  minimum  resolvable  wave  angle  is 
dependent  on  cos-1 the  possible  maximum  resolvable  wave  angle  is  dependent 
on  the  sample  interval,  which  will  be  discussed  in  the  next  section. 

The  above  derivation  is  based  on  the  assumption  of  the  time  sequential  wave 
cut  data.  In  some  applications,  however,  the  path  of  a  wave  cut  may  be  described 
by  a  path  equation  y  =  yo  +  tanas,  for  instance,  the  one-dimensional  wave  data 
is  cut  from  two-dimensional  wave  data  in  the  x  and  y  spatial  domain.  For  this 
case,  the  wave  elevation  is  given  in  the  form  of  fj(x)  =  rj(x,y(x)),  its  spectrum 
can  be  obtained  by  talcing  the  Fourier  transform  with  respect  to  x.  Then  similar 
procedures  to  calculate  the  wave  amplitude  function  can  proceed.  For  instance,  the 
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wave  amplitude  function  for  a  wave  cut  rj(x,  y0)  is  given  by 

Mff)  =  ~  (3.50) 


U  — 


=  0 


(3.51) 


2ir  U2  cos  9 

where  yo  is  a  constant  and  f/(u)  is  Fourier  transform  with  respect  to  variable  x. 

Similar  to  the  two-dimensional  case,  it  is  assumed  in  the  above  discussion  that  the 
truncation  error  caused  by  finite  data  length  can  be  neglected.  In  some  situations, 
however,  the  truncation  error  is  too  large  to  be  neglected,  and  thus  some  remedy 
methods  are  needed.  One  method  is  to  extend  the  truncated  wave  cut  according  to 
the  theoretical  ship  wave  asymptotic  behavior  as  given  in  [20].  Another  possibility, 
which  will  be  discussed  here,  is  the  use  of  an  inversion  technique  to  recover  the 
wave  amplitude  function  from  truncated  wave  spectra,  which  is  similar  to  the  two- 
dimensional  case  discussed  in  Section  3.1. 

In  the  second  method,  it  is  assumed  that  the  data  length  is  T,  and  that  the 
truncated  wave  cut  is  represented  by  fjr(t)  =  rj(t)  gr{t  —  f),  where  the  gate  function 
is  defined  as  gr(t)  =  1  while  ^  <  t  <  j  and  zero  otherwise.  The  Fourier  transform 
of  the  one-dimensional  gate  function  is  T  Sa(vT f) .  By  talcing  the  Fourier  transform 
of  t}t,  the  spectrum  of  the  truncated  wave  cut  is  given  in  the  form 

HT{f)  =  Te~W  y*  }*♦  eW*WSa{*T(f  -  |$(0)|))  dO  (3.52) 


When  Hr(f)  is  known,  A(9)  cam  be  found  by  solving  this  integral  equation.  Though, 
the  cadculation  will  be  more  complicated  than  the  method  neglecting  truncated  er¬ 
rors. 

So  far,  the  first  section  and  this  section  have  discussed  the  methods  to  recover 
the  wave  amplitude  function  from  either  one-  or  two-dimensional  wave  spectra.  Gen¬ 
erally  speaking,  the  calculation  for  one-dimensional  data  is  simpler  and  performance 
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is  better  than  those  for  two-dimensional  data  if  the  one-dimensional  data  contain 
enough  required  ship  information.  The  reason  that  the  accuracy  may  be  degraded 
in  two-dimensional  discrete  cases  is  that  the  spectrum  consists  of  discrete  sample 
pixels,  the  impulse  on  the  locus  can  not  be  always  located  on  these  sampling  points, 
and  a  slight  deviation  from  the  locus  may  result  in  a  large  error  for  that  impulse 
value.  However,  the  advantage  of  using  two-dimensional  data  is  that  it  is  much  eas¬ 
ier  to  remove  the  background  noise,  such  as  a  rough  wind  generated  wave,  and  it 
is  also  easier  to  estimate  a  ship’s  speed  and  direction  from  two-dimensional  wave 
fields  than  from  one-dimensional  wave  cuts.  For  the  above  reasons,  it  is  suggested 
that  the  signal  processing  and  the  estimation  of  ship  speed  and  direction  proceed 
in  the  two-dimensional  basis,  but  the  wave  amplitude  function  be  recovered  from 
one-dimensional  data  that  are  extracted  from  the  two-dimensional  data. 

3.3  Effect  of  Data  Sampling  Intervals 

The  wave  amplitude  function  is  a  function  of  the  wave  angle.  Theoretically,  the 
wave  angle  ranges  from  — *  to  *  for  ship  generated  surface  waves.  In  most  situations, 
however,  the  data  we  obtained  are  discrete,  and  the  range  of  the  wave  angle  is 
dependent  on  the  sample  interval  and  the  ship  speed  as  will  be  seen  below.  The 
following  two  subsections  discuss  the  effect  of  sampling  intervals  on  the  maximum 
resolvable  wave  angles  in  one  and  two-dimensional  cases,  based  on  Nyquist’s  sampling 
theory.  The  final  subsection  discusses  the  relationship  between  the  wave  sampling 
intervals  for  a  reed  ship  and  its  model,  because  ship  models  and  tow  tank  experimental 
data  are  usually  used. 
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3.3.1  Maximum  Resolvable  Wave  Angles  in  1-D  Cases 

The  maximum  resolvable  wave  angle  from  a  one-dimensional  wave  spectrum  is 
determined  in  this  subsection,  based  on  Nyquist’s  sampling  theory  and  the  relation 
between  the  wave  angle  and  wavelength. 

To  apply  Nyquist’s  sampling  theory  to  the  ship  wave  sampling  problem,  the 
wavelength  expression  is  given  first  here.  For  a  ship  moving  at  speed  U,  the  wave 
number  K  and  wave  angle  6  have  the  following  relation  from  Chapter  2: 

K  —  ^ 

U2  cos2#' 

According  to  the  definition  of  wavelength  A  and  (3.53),  it  follows  that 

.  A  2lf 

X  ~  K 

=  —u2cos2e 

9 

=  0.641I/2  cos20  .  (3.54) 

According  to  Nyquist’s  sampling  theory,  in  order  to  reconstruct  a  signal  from  its 
sampling  values  without  aliasing  error,  the  sampling  interval  A  in  the  spatial  domain 
must  be  such  that 

A  <  ^Amin  (3.55) 

where  Am,n  is  the  shortest  wavelength  that  the  signal  contains.  If  the  data  sampling 
is  taken  in  the  time  domain  mentioned  in  Section  3.2,  then  A  =  UpAt.  Applying 
this  sampling  criterion  to  the  above  ship  wave  problem,  it  follows  that 

A  <  —U2  cos2  6  .  (3.56) 

9 

Some  comments  can  be  made  here  in  terms  of  (3.56).  Since  the  wavelength  ranges 
from  0  to  22 (72,  corresponding  to  the  wave  angle  in  [— f,§],  the  sampling  interval 
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must  be  from  0  to  *U2.  Thus,  the  sampling  interval  must  approach  zero,  in  order  to 
cover  the  wave  angle  approaching  ±-|.  In  practice,  the  determination  of  the  sampling 
interval  depend*  on  many  other  facts,  and  the  interval  can  not  be  very  small.  For 
any  sampling  interval  A  larger  than  zero,  the  signal  components  with  the  wave  angle 
close  to  will  inevitably  have  some  distortion.  However,  if  the  intensities  of  these 
components  are  very  small  when  the  wave  angle  approaches  this  distortion  may 
be  neglected  in  real  applications. 

If  the  sampling  interval  A  is  given,  then  the  minimum  wavelength  which  can  be 
resolved  from  the  sampling  signal  is  determined  by  Amtn  =  2A;  additionally,  if  the 
ship  speed  U  is  also  given,  the  maximum  wave  angle  which  can  be  resolved  from  the 
sampling  signal  is  given  from  (3.56)  by 

)  (3-57) 

Thus,  the  available  signal  components  in  the  wave  spectrum  will  be  in  the  range  of 
0  =  0min  to  0  =  Omax,  where  0mi„  has  been  discussed  in  Section  3.2  and  is  given  by 
cos'*1  2 Ji^uup  if  A  =  UpAt  is  considered. 

For  easy  reference  with  different  ship  speed  parameters,  the  curve  of  maximum 
resolvable  wave  angles  from  one-dimensional  wave  spectra  versus  jfj  is  plotted  in 
Figure  3.5.  As  an  example  for  a  one-dimensional  case,  given  A  =  12.5  meters  and 
U  =  10  meters/second,  it  follows  the  maximum  resolvable  wave  angle  0max  =  51.36°, 
and  the  minimum  resolvable  wavelength  Amjn  =  25  meters. 

3.3.2  Maximum  Resolvable  Wave  Angles  in  2-D  Cases 

The  above  discussion  about  one-dimensional  cases  can  be  extended  directly  to 
two-dimensional  cases.  In  order  to  reconstruct  a  two-dimensional  signal  from  its 
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Figure  3.5:  Maximum  resolvable  wave  angle  from  1-D  or  2-D  wave  spectra  versus 

sampling  values  without  aliasing  error,  the  following  conditions  must  be  satisfied, 
according  to  Nyquist’s  sampling  theory: 


A*  S  2^*™"  — 


2u 


-  2Avm'"  “  2u, 


max 

l 


(3.58) 


where  A*  and  Ay  are  the  spatial  sampling  intervals  in  the  x—  and  y— directions, 
respectively.  AXmin  =  and  A,mjB  =  are  the  minimum  wavelengths  contained 
in  the  two-dimensional  signal  in  the  x—  and  y— directions.  umax  and  vmax  denote 
the  maximum  spatial  frequency  components  that  the  signal  contains.  If  sampling 
intervals  Ax  and  Ay  are  given,  the  spectrum  components  that  can  be  presented  are 
in  the  range  [— um,um]  and  [— vm,um]  in  the  spatial  frequency  (u,t?)  plane,  where 


Ur 


=:  4-  and  vm  =  At- 


2A, 


m  —  2A» 
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According  to  the  results  discussed  in  Section  3.1,  the  non-zero  spectrum  compo¬ 
nents  for  ship  waves  are  always  located  on  a  locus  on  the  (u,  v)  spatial  frequency 


plane,  as  represented  by 


2kU2  cos  0 
g  sin# 

2 irlP  cos2  0 


(3.59) 

(3.60) 


The  locus  will  have  an  intersection  point  with  the  edge  u  =  um  or  v  =  vm.  The 
maximum  resolvable  wave  angle  can  be  determined  from  the  intersection  point.  If 
the  locus  has  an  intersection  point  with  u  =  um,  then  the  maximum  resolvable  wave 
angle  can  be  obtained  by  solving  u  =  um  and  (3.59) . 


0  =  cos-1! - — - 1 

'  27TUmf/2 

=  ““(flf)  ■ 


(3.61) 


If  the  locus  has  a  intersection  point  with  v  =  vm,  then  the  maximum  resolvable  wave 
angle  can  be  obtained  by  solving  v  =  vm  and 


Sm  1  ^  2  ^  2irU*vm  +  \/(  2 irU3vm  )2  +  4  ^ 


(3.62) 


In  most  cases,  the  sampling  intervals  in  the  x—  and  y— directions  are  set  to  be  the 
same,  i.e.  A*  =  Ay  =  A.  It  can  be  proved  that  in  the  first  quadrant  of  the  (u,v) 
plane,  the  locus  is  under  the  line  u  =  t;forO<0<^,  and  is  above  the  line  u  =  v 
for  \  <  9  <  f.  From  (3.59)  and  (3.60), 


thus, 


U-“=2 


if  0  <  0  <  - 
4 


u  >  v 


(3.63) 
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it  <  v  if  ^  <  0  <  ^  .  (3.64) 

From  (3.64),  it  can  be  concluded  for  the  cases  of  A*  =  Av  =  A  that  when 
9 max  <  the  locus  has  an  intersection  point  with  u  =  uTO  and  when  \  <  6 max  <  f , 
the  locus  has  an  intersection  point  with  t;  =  vm,  and  A  can  be  expressed  from  (3.60) 
in  the  form 

jU2  COS  9 max  if  0  <  9max  <  f 

A  =  < 

jU2  COS  9max  COt  9 max  if  f  <  9max  <  f 

(3.65) 

Figure  3.5  also  gives  the  maximum  resolvable  wave  angle  curves  versus  the  ratio  of 
the  ship  speed  U  and  the  sampling  interval  A  =  A*  =  Ay.  From  Figure  3.5,  we  find 
that  the  sampling  intervals  must  be  small  enough  to  recover  the  wave  components 
with  large  wave  angles.  Considering  the  above  example  again  for  a  two-dimensional 
case,  given  A  =  12.5  meters,  U  —  10  meters/second,  then  9max  =  55.47°  in  terms 
formula  (3.62). 


3.3.3  Sampling  Intervals  for  a  Real  Ship  and  its  Model 


This  subsection  discusses  the  relation  between  the  wave  sampling  intervals  for  a 
real  ship  and  its  model.  For  this  purpose,  a  non-dimensional  parameter  Fn,  called 
the  Froude  number,  is  introduced  here.  With  physical  length  L,  speed  U  and  gravi¬ 
tational  acceleration  g,  Fn  is  defined  by 


(3.66) 


It  is  widely  used  in  the  ship  wave  resistance  analysis.  According  to  studies  of  ship 
wave  resistance,  the  wave  resistance  of  two  geosims  with  the  same  hull  shape  are  the 
same  when  their  Froude  numbers  are  equal  [21]. 
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If  applying  the  above  result  to  a  real  ship  with  length  L  and  moving  speed  U  and 
its  model  with  the  same  hull  shape  but  length  Lm  and  moving  speed  Um,  then  the 
following  ratio  can  be  obtained  when  their  Froude  numbers  are  equal: 

V 1  L 

To  express  the  sampling  interval  relation,  let  the  sampling  intervals  corresponding 
to  the  real  ship  and  its  model  be  denoted  by  A  and  Am,  respectively.  The  ratio  of 
the  two  sampling  intervals  can  be  connected  to  the  above  ratio  in  the  relation 

A m  =  Ul=Ln 

A  U*  L 

for  either  the  one-dimensional  case  in  (3.57)  or  the  two-dimensional  case  in  (3.65). 
This  tells  us  that  when  the  Froude  number  is  kept  the  same,  the  sampling  interval 
corresponding  to  the  ship  model  can  be  taken  as  Am  =  ^A. 

So  far,  this  section  has  discussed  the  effect  of  sampling  intervals  on  the  maximum 
resolvable  wave  angles  in  one  and  two-dimensional  cases,  and  the  relation  between 
the  wave  sampling  intervals  for  a  real  ship  and  its  model.  In  summary,  from  the 
viewpoint  of  analysis  of  ship  wave  spectra  and  extraction  of  ship  geometry  informa¬ 
tion,  the  determination  of  sampling  intervals  depends  mainly  on  the  ship  speed  and 
ship  length.  The  smaller  the  sampling  interval,  the  smaller  the  distortion  of  wave 
spectra,  the  larger  the  maximum  resolvable  wave  angle.  Additionally,  it  will  be  seen 
in  Chapter  4  that  the  smaller  the  sampling  interval,  the  more  periodic  zero  points  are 
available  in  the  wave  amplitude  function.  Therefore,  the  sampling  interval  should  be 
selected  as  small  as  possible.  In  practice,  however,  the  resolution  and  properties  of 
data  acquisition  systems  and  their  operation  position,  such  as  SAR,  greatly  limit  the 
small  sampling  intervals  to  be  used.  Other  limitations  may  be  the  data  storage  and 
data  processing,  but  they  are  less  critical  compared  to  the  former.  For  the  wave  gen- 
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erated  by  a  ship  model,  the  sampling  interval  can  be  taken  as  the  one  proportional 
to  the  ratio  of  the  lengths  of  the  model  and  real  ship. 


CHAPTER  IV 


ESTIMATION  OF  A  MOVING  SHIP’S  SPEED 

AND  DIRECTION 


From  the  viewpoint  of  remotely  sensing  a  ship  moving  in  the  open  ocean,  inter¬ 
esting  problems  exist  related  to  the  detection  of  the  ship’s  actual  presence,  and  the 
acquisition  of  its  dynamic  and  static  information,  for  instance,  the  ship’s  direction 
and  speed  and  the  ship’s  size  and  hull  shape.  These  problems  will  be  discussed  in 
the  following  chapters  based  on  the  knowledge  given  in  the  previous  chapters.  This 
chapter  focuses  on  the  estimation  of  a  ship’s  direction  and  speed  from  one  dimen¬ 
sional  and  two  dimensional  wave  spectra.  Before  this  discussion,  the  problem  of  the 
presence  of  a  moving  ship  in  ambient  ocean  waves  is  briefly  discussed. 

4.1  Presence  of  a  Moving  Ship  in  Ambient  Ocean  Waves 

Although  a  ship’s  length  is  bounded  within  a  range  of  values,  the  wake  it  generates 
in  the  open  ocean  may  extend  for  many  tens  of  kilometers.  In  the  indirect  methods, 
the  detection  of  a  ship  and  its  related  characteristics  is  obtained  by  measuring  ship 
generated  waves  or  their  spectra. 

One  important  feature  of  the  ship  wake,  different  from  that  of  ambient  ocean 
waves,  is  its  wave  spectrum.  As  analyzed  in  Chapter  2  the  spectrum  of  the  complex¬ 
valued  wave  elevation,  H( u,u),  has  one  locus  on  the  spectrum  diagram,  and  the 
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Figure  4.2:  (a)  Pseudo  image  of  the  Quapaw  ship  wave  elevation  in  a  random  sine 
ambient  wave,  (b)  Pseudo  image  of  the  Fourier  transform  of  the  wave 
elevation  from  (a). 
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spectrum  of  the  real-valued  wave  elevation,  H(u,v ),  has  two  loci  on  the  spectrum 
diagram.  Because  of  the  distinct  spectrum  characteristics,  it  is  usually  easier  to 
detect  the  presence  of  a  moving  ship  from  the  sea  noise  background  in  the  spectral 
domain  them  in  the  spatial  domain. 

To  understand  this  point,  let  us  examine  a  simple  example  of  a  ship  wave  plus  a 
random  sine  ambient  wave.  This  is  meant  to  simulate  a  ship  in  a  swell  background. 
The  ship  wave  without  any  ambient  waves  and  background  noise  is  shown  as  a 
pseudo  image  in  Figure  4.1(a)  together  with  its  FFT  spectrum  in  Figure  4.1(b). 
In  Figure  4.1(a),  the  origin  of  the  coordinate  system  has  a  translation  and  rotation 
relative  to  the  ship  center,  the  origin  of  the  ship  reference  system  defined  in  Chapter  2. 
This  difference  results  in  a  rotation  of  the  loci  on  the  spectrum  diagram,  but  it  dose 
not  change  the  shape  of  the  loci.  Further  discussion  about  it  will  be  given  in  the 
next  section. 

This  ship  wake  is  calculated  using  WAVEAMP,  a  program  to  compute  the  Kelvin 
wave  elevation  [29],  for  a  1:12  scale  model  of  a  seagoing  tug,  the  USS  Quapaw,  which 
has  a  length  of  4.953  meters  and  a  speed  U  of  2.229  meters/second  with  an  angle  of 
10°  relative  to  the  x-axis.  The  ship  wave  height  has  a  maximum  value  0.231  ,  mean 
0.007  and  standard  deviation  0.013  meters.  The  ship  wave  involved  in  a  random  sine 
ambient  wave  is  displayed  in  Figure  4.2(a)  together  with  its  FFT  spectrum  in  Figure 
4.2(b).  The  random  sine  wave  has  a  simple  model,  Ai,sin(Kxx  +  Kyy),  where  Ab  Kx 
and  Ky  are  random  variables  generated  point  by  point  by  a  computer  program.  Ab 
originally  has  a  Gaussian  distribution  with  mean  0.05  meters  and  standard  deviation 
0.05  meters,  denoted  as  A/*(0.05, 0.05),  and  Kx  and  Kv  originally  have  a  Gaussian 
distribution  .V(3.14, 1.0)  in  rad./meter.  They  are  smoothed  using  a  median  filter  with 
a  9-point  window  size,  equivalent  to  1.8  by  1.8  meters.  The  smoothed  Ab,  Kx  and 
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Ky  have  means  close  to  their  original  means,  but  have  different  standard  deviations 
of  0.014  meters,  0.28  rad. /me ter  and  0.28  rad./meter,  respectively.  Finally,  they  are 
used  to  calculate  the  random  sine  wave  and  are  added  to  the  ship  wave. 

From  this  example,  it  can  be  found  that  in  the  spatial  domain,  the  ship  wave, 
particularly  the  wave  on  the  left  of  the  ship,  has  been  corrupted  by  the  random  sine 
wave  because  of  their  close  wave  direction,  but  the  loci  can  be  still  recognized  clearly 
from  the  spectrum.  In  real  situations  with  severe  background  noise,  conventional  or 
special  signal  processing  may  be  used  to  enhance  the  desired  ship  wave  signal. 

4.2  Estimation  of  a  Ship’s  Speed  and  Direction  from  2-D 
Wave  Spectra 

This  section  discusses  the  estimation  of  a  moving  ship’s  speed  and  direction  from 
its  two  dimensional  wave  spectrum.  The  discussion  will  begin  with  two  kinds  of 
spatial  coordinate  systems  and  their  corresponding  spectrum  coordinate  systems, 
and  then  the  formulas  for  estimating  the  speed  and  direction  are  derived. 

In  the  following,  the  discussion  will  focus  on  the  estimation  from  the  magnitude 
of  a  Fourier  spectrum,  instead  of  the  one  from  a  power  spectrum,  since  a  power 
spectrum  and  the  magnitude  of  a  Fourier  spectrum  have  a  direct  relation  and  are 
equivalent  when  the  spectrum  locus  position  is  used  to  estimate  a  ship’s  speed  and 
direction. 

It  has  been  shown  in  Chapter  2  that  under  the  steady  state  assumption,  the  ship 
speed  U,  wave  angle  0  and  wave  number  K{9)  have  a  direct  relation 

■  <41> 

This  relation  indicates  that  the  ship  speed  depends  only  on  the  wave  number,  or 
the  wave  length  at  a  given  wave  angle.  Theoretically,  once  the  wave  number  K(9) 
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is  available,  the  determination  of  ship  speed  becomes  a  trivial  problem.  However, 
the  determination  of  the  wave  number  and  wave  angle  needs  prior  information  about 
the  ship’s  direction.  If  this  information  is  not  available,  then  the  problem  becomes 
complicated,  and  both  the  ship’s  direction  and  speed  need  to  be  determined  simul¬ 
taneously. 

To  determine  the  ship’s  direction,  two  coordinate  systems  in  the  spatial  domain, 
shown  in  Figure  4.3,  are  considered  in  this  section.  One  is  the  reference  system 
moving  with  a  ship  as  defined  in  Chapter  2;  another  is  the  image  coordinate  system 
whose  origin  is  the  imaging  center  and  the  positive  x-direction  is  the  sensor’s  direc¬ 
tion.  If  xoy  denotes  the  ship  reference  coordinate  system  and  xmomym  the  image 
coordinate  system,  then  their  relation  is  given  as 

xm  =  xm0  +  x  cos  a  —  y  sin  a  (4.2) 

ym  =  J/mo  +  x  sin  or  -f-  y  cos  a  (4.3) 

or 

x  =  (im  -  *mo)  COS  a  +  (ym  -  ym0 )  sin  a  (4.4) 

y  =  -(im  -  *mo)  sin  a  +  (ym  -  ymo)  cos  a  (4.5) 

where  (xmo,  ym o)  is  the  coordinate  of  the  origin  o  of  the  ship  reference  system  in  the 
image  coordinate  system  and  where  a  is  the  angle  between  axes  ox  and  omxm,  which 
represents  the  ship’s  direction  relative  to  the  sensor’s  direction.  With  the  relations, 
the  ship  wave  ?y(x,y)  expressed  in  the  ship  reference  system  can  be  expressed  in  the 
image  coordinate  system  as 

r)m{xm,ym)  ~  ?/[  (xm  -  xm0)  cos  a  +  (ym  -  ym0)  sin  a, 


-(Xm  -  xmo)sina  +  (ym  -  ym0)cosa  . 


(4.6) 


Figure  4.3:  Two  coordinate  systems  in  spatial  domain:  the  ship  reference  coordinate 
system  xoy  and  image  coordinate  system  xmomym 


Figure  4.4:  Two  coordinate  systems  in  frequency  domain:  the  ship  reference  coordi¬ 
nate  system  uov  and  image  coordinate  system  umomvm 
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Corresponding  to  the  the  above  spatial  domain  coordinate  relations,  the  spec¬ 
trum  domain  coordinate  relations  can  be  established  through  the  Fourier  transform 
relation.  The  Fourier  transform  of  ij(x,y)  has  been  given  in  (2.20),  i.e.,  H{u,v)  = 
T{  T)m(x,y)  }.  The  Fourier  transform  of  ym)  is  given  as 


Hm(um,vm)  =  ^F{Vm(xm,ym)}  =  [  [  Vm(xm,  ym)e~i2ir(UmXm+Vmym)  dxmdy> 

J —oo  J —oo 

/oo  too 

/  l[(xm  -  Xm0)  cos  a  +  (ym  -ym0)sina, 

-oo  J— oo 

-(xm  -  xm0)  sin  a  +  (ym  -  ym0)  cos  a  ]e~i2^UmXm+VmVm)dxmdym 
=  [  /°°  /°°  Tf(x,y)e-j2*[x(UmC~a+Vm'iaa)+y(‘-Umlina+VmCOaa)]dxdy] 

J —oo  J —OO 


-oo  j—  oo 

,«”j  2ir(ttmxmo+Vm  VmO  ) 


=  H(  um  cos  a  +  vm  sin  a,  — ttm  sin  a  -f  vmcosa  )e  i2ir(u™I"»o+,'’»»mo)  (4.7) 

From  the  above  relation  of  the  Fourier  transforms  Hm(um,vm)  and  H(u,v ),  the 
coordinate  relation  in  the  spectrum  domain  is  given  as 


u  =  um  cos  a  +  vm  sin  a  (4.8) 

v  =  —  vm  sin  a  +  vm  cos  a  (4.9) 


or 


um  =  ucosa  —  t/sina  (4-10) 

vm  =  u  sin  a  +  v  cos  a  .  (4.11) 


In  terms  of  the  above  expressions,  the  following  comments  can  be  made: 

(1)  The  spectrum  domain  coordinate  system  uov,  corresponding  to  xoy,  has  a  ro¬ 
tation,  with  an  angle  a,  relative  to  the  spectrum  domain  coordinate  system  umomvm, 
corresponding  to  xmomym.  Thus,  if  a  is  determined,  then  the  ship’s  direction  relative 
to  the  sensor’s  direction  can  be  obtained. 
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(2)  The  translation  (xmo,  ymo)  between  the  two  spatial  coordinate  systems  reflects 

only  on  the  phase  of  vm).  Thus,  it  is  not  important  to  know  the  translation 

since  the  wave  number  and  wave  angle  are  determined  from  the  magnitude  of  the 
Fourier  spectrum  or  the  power  spectrum. 

(3)  When  a  =  0,  then  | Hm(um,  t?m)|  =  vm)|,  the  wave  number  components 

and  the  spatial  frequencies  have  simple  relations,  i.e.,  Kx  =  2irum  and  Ky  =  2irvm, 
and  the  wave  angle  is  given  as  0  =  For  this  case,  the  speed  can  be  calculated 

by  directly  measuring  the  position  of  the  locus  points  in  the  spectrum  diagram. 

(4)  According  to  the  spectrum  coordinate  relation,  it  follows  txjj,  +  =  u2  4-  v2; 

thus,  the  wave  number  K  is  invariant  with  the  coordinate  system  transform,  that  is, 

K  —  2ir\Ju2n  +  =  2 vy/ti2  +  v2  .  (4-12) 


With  the  above  relations  and  the  conclusions,  the  general  formulas  to  determine 
a  moving  ship’s  direction  and  speed  are  derived  in  the  following.  First,  consider  the 
ship’s  wave  spectrum  locus  in  the  uov  coordinate  system.  From  (4.1),  Kcos29  =  fj; 
thus,  for  any  two  locus  points,  it  follows 


(4.13) 


where  the  subscript,  1  or  2,  indicates  that  the  wave  number  and  angle  are  obtained 
from  the  given  point  1  or  2.  Since  cos  9  =  y  =  (4.13)  can  be  rewritten  as 

^ui  =  .  (4.14) 

In  the  umomvm  coordinate  system,  (4.14)  becomes 

y/K^(umtcosa +  vmisina)  -  yfK[{um2  cos  a  +  vm2s\n  a)  (4.15) 


and,  thus,  the  angle  a  can  be  estimated  from 

.  .  _i  y/J^ium2  —  y/Klumi 

a~  VJr,vmt  -  ^Jr2vm, 


(4.16) 
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Similarly,  the  ship  speed  can  be  calculated  from 


9 

KjCos29j 


fT 


V  2ir  \umjCosa  +  vmj  sin  d| 


(4.17) 


where  j  =  1  or  2. 

In  real  situations,  the  ship  wave  spectrum  is  discrete,  and  thus  the  locus  points 
on  the  spectrum  diagram  will  not  always  exactly  locate  on  the  sampling  grids.  This 
results  in  errors  on  some  locus  points  when  they  are  read  from  the  spectrum  diagram. 
To  remove  this  effect  on  the  calculation  of  the  ship’s  speed  and  direction,  many  pairs 
of  locus  points  can  be  used  to  calculate  the  ship’s  direction  and  speed,  and  then  their 
average  is  taken  as  the  estimate  of  the  direction  and  speed.  Specifically,  consider 
there  are  M  pairs  of  locus  points  available.  The  angle  a  is  calculated  with  each  pair 
of  points,  and  the  average  of  the  calculated  angles  a*,  i  =  1, ...,  M  is  then  considered 
as  the  estimate  of  a,  i.e., 

1  M 

“=  M  52  •  (4-!8) 

The  average  of  f/,,  t  =  1,  ...,2 M  is  considered  to  be  the  estimate  of  the  ship  speed  , 
that  is, 


U  = 


2  Mi 


2  M 


«=i 


J_v  rr  vk 

2M  JrJ  V  27r  | umiCosa  +  vmi  sin  q| 


(4.19) 


The  formulas  for  estimating  a  moving  ship’s  direction  and  speed  from  a  two 
dimensional  wave  spectrum  have  been  derived  above.  The  scheme  for  the  estimation 
is  now  shown  in  Figure  4.5.  In  real  situations,  ship  waves  are  involved  in  a  random 
sea  background.  Thus,  a  wave  spectrum  contains  not  only  the  ship  wave  components 
but  also  the  noise  components.  To  remove  the  random  noise  and  other  undesired 


49 


Figure  4.5:  Scheme  for  estimating  a  moving  ship’s  direction  and  speed  from  its  wave 
spectrum 

components,  digital  processing  techniques  may  be  used.  The  clipper  shown  in  the 
figure  is  used  for  the  this  purpose. 

Theoretically,  the  two  dimensional  Fourier  transform  of  the  ship  generated  wave 
is  composed  of  6-functions;  thus,  there  are  many  infinite-size  impulses  located  on  the 
loci  of  the  spectrum  diagram  as  discussed  in  Chapter  3.  When  the  finite  wave  patch 
is  sensed  as  in  real  situations,  they  are  on  the  order  of  0(LC)  [6],  where  Lc  denotes 
the  characteristic  length  of  the  finite  patch.  For  the  case  of  high  signal  to  noise 
ratio,  a  simple  processing  method  can  be  used.  For  instance,  a  clipper  is  used  to 
remove  the  background  noise  components.  This  processing  is  helpful  in  determining 
the  position  of  each  locus  point.  From  the  positions  of  locus  points,  the 

ship’s  direction  and  speed  finally  are  calculated.  The  algorithm  can  be  implemented 
in  software  with  a  fast  and  accurate  estimation  performance. 

To  demonstrate  the  above  method,  consider  here  an  example  of  the  Quapaw’s 
wave  elevation  field,  shown  in  Figure  4.1(a).  It  is  assumed  in  the  calculation  that 
the  ship’s  direction  and  speed  are  a  =  10°  and  U  =  2.229  meters/second,  and  that 
the  sampling  intervals  in  the  x-  and  y-directions  are  0.2  meters.  The  ship’s  direction 
and  speed  are  estimated  from  the  spectrum  of  the  wave  field,  whose  contour  plot  is 
shown  in  Figure  4.6.  Note  that  the  subscript  “m”  has  been  used  in  this  figure  to 
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u»(1M) 

Figure  4.6:  Contour  plot  of  the  Fourier  transform  vm)  of  the  Quapaw’s  wave 

elevation. 

emphasize  the  image  coordinate  system.  A  total  of  183  locus  points  are  evaluated 
for  the  right  locus  on  the  spectrum  with  a  threshold  of  /xjy  +  <th,  where  hh  and  an 
are  the  mean  and  standard  deviation  of  the  spectrum  intensity.  According  to  the 
computer  calculation,  the  estimated  direction  is  a  =  9.999°  with  a  relative  error  of 
0.006%  and  a  r.m.s.  error  of  0.233°;  the  estimated  speed  is  U  =  2.230  meters/second 
with  a  relative  error  of  0.03%  and  a  r.m.s.  error  of  0.019  meters/second. 

4.3  Estimation  of  a  Ship’s  Speed  and  Direction  from  1-D 
Wave  Spectra 

This  section  discusses  the  estimation  of  a  moving  ship’s  speed  and  direction  from 
its  one  dimensional  wave  spectrum.  The  expression  for  the  spectrum  of  a  wave  cut 
making  an  angle  with  the  positive  x-axis  has  been  given  in  (3.40).  The  wave  cut 
spectrum  has  two  peaks  under  certain  conditions,  and  the  frequency  positions  of  the 
two  peaks  can  be  used  to  estimate  the  ship’s  direction  and  speed,  as  suggested  by 
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Tuck  et  al  [6].  For  this  one-cut  method,  the  formulas  for  estimating  the  direction 
and  speed  are  given  here  first,  and  then  another  method,  called  the  two-cut  method, 
is  introduced  for  estimating  the  ship  speed  from  two  wave  cuts  parallel  to  the  ship’s 
central  line.  Examples  are  also  given  to  demonstrate  the  two  methods. 

4.3.1  One-cut  Method 

The  discussion  begins  with  the  wave  cut  spectrum  and  its  properties.  The  spec¬ 
trum  expression  for  a  wave  cut  have  been  given  in  (3.40),  that  is, 

with  a  constraint  condition  for  0,  /  —  |$(0)|  =  0,  as  in  (3.41).  This  condition 
indicates  9  is  a  function  of  frequency  /;  thus,  there  may  exist  some  frequencies  such 
that  $  (0)  =  0,  and  hence  there  may  exist  some  singularities  for  the  spectrum.  For 
the  finite  length  wave  cut,  this  will  cause  some  sharp  peaks  on  its  spectrum  diagram. 
The  peak  height  is  proportional  to  the  square  root  of  the  data  record  length  [6].  The 
frequency  points  of  the  peaks  on  the  wave  cut  spectrum  can  be  determined  by  two 
equations,  (3.41)  and 

$'(0)  =  0  .  (4.20) 

Here,  $  depends  on  the  wave  angle  9  as  well  as  the  wave  cut  angle  a.  Substituting 
the  approximation  of  $  given  in  (3.43),  $(0)  «  into  (4.20)  and  then 

solving  the  resulting  equation  together  with  /  — |$(0)|  =  0  yields  the  relation  between 
9  and  a: 

sin  29 


a  =  tan 


cos  20  —  3 


(4.21) 
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The  curve  of  a  versus  9  is  plotted  in  Figure  4.7  (a).  It  is  found  from  the  figure  or 
calculation  that  |s:n-1||(=  19.5°)  is  the  maximum  wave  cut  angle  for  the  peaks  to 
exist. 


Figure  4.7:  Relations  between  the  wave  angle  9,  cut  angle  a  and  frequency  fv  at 
peaks,  (a)  a  versus  9 ;  (b)  versus  0;  (b)  versus  a;  (b)  ^ 

versus  a. 


With  the  above  relation  of  9  and  a,  the  special  frequency  points  fv  are  obtained 
from  /  —  |$(0) |  =  0  in  the  form 

2*/P  _  1 


KqUp  cos  9\J  4  —  3  cos2  6 


(4.22) 


The  frequencies  also  can  be  directly  expressed  in  a  function  of  the  wave  cut  angle  a: 

2* fP  _ 


KoUp  \\  4  - 


(4.23) 


where  Ko  =  and  Up  and  U  are  the  speeds  of  the  sensor  and  ship,  respectively, 
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and  where 


n,2(a)  = 


(1  —  8  tan2  a)  ±  3\/l  —  8  tan2  a 


(4.24) 


2(1  +  tan2  a) 

represents  the  value  rj  with  or  r2  with  ”  before  the  sign  of  squaxe  root. 
Since  ri  and  r2  are  real,  1  —  8 tan1  a  >  0,  which  is  equivalent  to  a  <  |sm-1||  as 
asserted  above.  In  terms  of  (4.23)  and  (4.24),  there  axe  one  or  two  frequency  points, 
corresponding  to  n  and  r2,  where  the  peaks  appear.  When  a  =  0,  there  is  a  finite 
frequency  point  fpl  and  an  infinite  frequency  point  /p2  =  oo;  when  0  <  a  <  |sm-1||, 
there  are  two  finite  frequency  points,  fpl  and  /p2;  when  a  =  |sin-1i|  there  is  only 
one  finite  frequency  point  fp\  =  /p2. 

From  the  above  special  frequency  points,  the  ship’s  direction  is  estimated  first 
and  then  the  speed  is  calculated.  When  0  <  a  <  |sm-1l|,  the  direction  is  estimated 
by  calculating  the  ratio  of  the  above  special  frequencies,  i.e., 


fpi 

fp2 


4  -  r*(or) 


(4.25) 


\J  4  —  r2(a 

The  frequencies  /pl  and  fp2  are  found  from  the  wave  cut  spectrum,  then  the  wave 
cut  angle  is  obtained  by  solving  (4.25),  and  finally  the  ship  speed  is  calculated  from 
(4.23)  by  noting  the  relation  of  Kq  with  the  speed,  i.e.,  Kq  =  Figure  4.7  (d) 
shows  the  curve  of  the  ratio  ^  versus  the  wave  cut  angle  a.  When  a  =  0,  the  ship 
speed  is  directly  estimated  from  /pi  with  the  formula 

2*7Pi 


K0U, 


=  1 


(4.26) 


To  demonstrate  the  above  one-cut  method,  consider  here  an  example  in  which 
the  wave  elevation  cut,  shown  in  Figure  4.8(a)  is  calculated  using  WAVEAMP.  The 
ship  model  is  the  same  as  in  the  above  two  dimensional  case.  Each  cut  has  256  data 
points  and  a  sampling  interval  of  0.001  seconds.  It  is  assumed  that  the  ship  has  a 


IFFT  of  wave  elevation  I 
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Figure  4.8:  (a)  Wave  elevation  cut  calculated  from  WAVEAMP  for  the  Quapaw 
model  with  wave  cut  angle  a  =  10°,  speed  U  =  2.229  m/s,  and 
(®mOiJ/mo)  =  (-48.5087,-20.5690);  (b)  Magnitude  of  the  Fourier  trans¬ 
form  (dotted  line)  of  the  wave  elevation  from  (a),  and  two  peaks  (solid 
lines)  where  /  —  |$(0)|  =  0. 


55 


direction  of  a  =  10°  relative  to  the  sensor’s  direction,  and  that  the  sensor’s  speed 
is  Up  =  199.219  meters/second.  The  frequency  points  at  peaks  are  obtained  from 
the  spectrum,  which  is  the  dotted  line  in  Figure  4.8(b).  The  two  desired  peaks  have 
been  detected  and  are  shown  in  solid  lines  in  the  figure.  Since  the  two  peaks  are 
located  close  together,  the  peak  P2  has  a  little  left  shift  in  the  spectrum  shown  in  the 
dotted  line.  Thus,  the  peak  P3  has  to  be  determined  after  removing  the  peak  Pi  from 
the  spectrum  using  a  bandpass  filter  with  good  selectivity.  The  filtering,  position 
detection  and  all  other  calculation  can  be  automatically  completed  using  a  computer 
program.  According  to  the  computer  calculation,  the  detected  position  of  the  two 
peaks  are  fp\  =  62.5  Hx  and  fpi  =  95.7  Hz\  the  estimated  direction  is  ot  =  10.019° 
with  a  relative  error  of  0.19%;  the  estimated  speed  is  U  =  2.195  meters/second  with 
a  relative  error  of  1.51%. 


4.3.2  Two-cut  Method 

When  a  wave  cut  is  parallel  to  the  ship’s  central  line,  the  wave  cut  angle  is  equal 
to  zero.  For  this  case,  an  alternative  method  can  be  used,  which  is  called  here  the 
two-cut  method  because  two  wave  cuts  parallel  to  the  ship’s  central  lines  are  used.  In 
this  method,  the  ship  speed  is  estimated  from  the  relative  phase  difference  between 
the  Fourier  transforms  of  the  two  wave  cuts.  The  derivation  and  example  are  now 
given  below. 

First,  recall  the  relation  between  the  wave  amplitude  function  and  the  Fourier 
transform  of  a  wave  cut  with  a  =  0,  which  are  given  in  (3.46)  and  (3.47): 


4(0)  _  ^(j)i^vW]ej(iCI(9)xo+A-v(e)I«) 
2ir 


(4.27) 


/- 


2vU  cos  6 


=  0  . 


(4.28) 
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Now,  consider  two  parallel  wave  cuts  with  cut  angle  a  =  0,  distance  Ay  =  ya  —  Vi 
and  the  same  starting  point  xo  in  the  x-direction.  The  Fourier  transforms  of  the  two 
wave  cuts  are  H\(f)  and  /^(Z).  Since  the  wave  amplitude  function  A{6)  is  supposed 
to  be  the  same  for  the  two  wave  cuts,  using  the  above  relation  and  dividing  H^if) 
by  H\(f)  leads  the  following  relation 


Ml = ei a® 

&(/) 


(4.29) 


where  the  relation  Ky  =  has  been  used.  If  the  phase  difference  of  /^(Z)  and 

Hi(f)  is  represented  by  A then  the  phase  relation  in  (4.29)  can  be  written  in  the 
form 


9  sin  0(f)  _  A <j>H 
U2  cos2  0(f)  fAy 

From  (4.28),  cos0(/)  =  -^Ju-  Substituting  this  relation  into  (4.30)  yields 

(2  «)4U2n  ,  9 

g2  1  k2x fUl  J  fAy  ' 


(4.30) 


(4-31) 


Note  that  from  the  above  expression,  approaches  1  as  A<f>n  approaches  0.  Thus, 
once  we  can  find  the  frequency  /m<„  corresponding  to  the  minimum  value  of  |A<£j/|, 
then  the  ship  speed  can  b.e  found  from 


2tT  Zm»'n 

To  demonstrate  this  two-cut  method,  consider  an  example  in  which  the  ship  wave 
elevation  cuts  were  measured  by  three  capacitance  wave  probes  when  the  Quapaw 
model  was  towed  in  a  tow  tank  l.  Three  wave  cuts  were  obtained  for  each  run. 
The  wave  elevation  of  two  runs,  RUN3  and  RUN5,  are  shown  in  Figure  4.9.  As 

an  example,  Figure  4.10  shows  the  magnitude  of  the  spectra  of  wave  cuts  RUN5-B 

^he  experiments  were  made  by  Ship  Hydrodynamics  Laboratory,  Department  of  Naval  Archi¬ 
tecture  and  Marine  Engineering,  the  University  of  Michigan  in  October,  1990. 


Figure  4.9:  Wave  elevation  cuts  from  tow  tank  experiments  for  the  Quapaw  model. 


I 
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f  (Hz) 


Figure  4.10:  Magnitude  of  the  FFT  of  wave  cuts,  RUN5-B  and  RUN5-C,  and  their 
phase  difference  A <f>n. 

and  RUN5-C  and  their  phase  difference.  Here,  A <j>H  has  been  processed  using  a 
median  filter  with  a  3-point  window  size  for  detecting  the  minimum  point  of  the 
phase  difference.  The  estimation  results  are  listed  in  Table  4.1. 


So  far,  the  methods  for  estimating  a  moving  ship’s  direction  and  speed  from  one 
and  two  dimensional  wave  spectra  have  been  discussed.  Comparing  these  methods, 
the  two  dimensional  method  has  three  primary  advantages.  First,  it  has  no  limitation 
on  the  ship’s  moving  direction,  except  for  the  180°  ambiguity  that  results  if  no  further 
prior  information  is  used.  This  is  opposed  to  the  one-cut  method,  which  is  suitable 
only  for  the  cut  angles,  a,  between  [—19.5°,  19.5°],  but  is  unable  to  tell  the  positive 


59 


RUNS 

RUN5 

cuts 

U  (m/s) 

error  (%) 

U  (m/s) 

error  (%) 

a  &  b 

2.283 

2.4 

2.283 

mm 

a  &  c 

2.283 

2.4 

2.283 

2.4 

b  &  c 

2.283 

mm 

2.179 

2.2 

Table  4.1:  Ship  speed  estimated  from  the  Quapaw’s  wave  cuts  using  the  two-cut 
method. 

or  negative  angle,  and  may  result  in  false  detection  if  obvious  peaks  also  exists  when 
|a|  >  19.5°.  Second,  the  two  dimensional  method  works  well  even  in  the  presence 
of  ambient  waves  and  background  noise  because  of  the  spectrum  feature  of  ship 
generated  waves.  In  the  one-cut  method,  the  ship  wave  signal  can  be  easily  corrupted 
by  ambient  waves  and  background  noise,  and  thus  it  may  result  in  false  detection 
of  the  peak  position.  Third,  the  two  dimensional  method  appears  to  achieve  more 
accurate  estimation  results  than  the  one  dimensional  methods.  Because  of  these 
reasons,  the  two  dimensional  method  should  be  always  considered  first  when  two 
dimensional  spectra  are  available.  It  has  recently  become  possible  to  obtain  these 
data  from  air-borne  or  space-borne  radar  systems  or  other  modem  remote  sensing 
techniques. 


CHAPTER  V 


ESTIMATION  OF  SHIP  LENGTH 

This  chapter  gives  a  detailed  discussion  of  the  estimation  of  a  ship’s  length  from 
its  amplitude  function.  The  ship  length  is  am  important  quantity  to  be  estimated 
in  remotely  sensing  ship  characteristics,  and  it  is  also  an  important  parameter  in 
the  further  estimation  of  ship  hull.  In  recovering  a  ship’s  hull  shape  from  the  wave 
amplitude  function,  an  inversion  problem  is  involved,  that  is,  an  integral  equation 
must  be  solved,  which  will  be  discussed  in  the  next  chapters.  The  integral  limits 
along  the  x-direction  are  specified  by  the  ship  length.  Therefore,  the  accuracy  of 
the  recovered  ship  hull  shape  will  depend,  to  a  great  extent,  on  the  accuracy  of  the 
estimated  ship  length. 

The  principle  of  the  estimation  of  ship  length  is  that  there  is  the  relation  between 
a  ship’s  hull  and  its  wave  amplitude  function,  and  that  ship  length  information  is 
contained  in  the  observable  periodic  character  of  the  wave  amplitude  function.  This 
character  can  be  found  not  only  in  the  read  and  imaginary  parts  of  the  wave  amplitude 
function,  but  also  in  the  magnitude  of  the  wave  amplitude  function. 

In  this  chapter,  a  theoretical  model  of  the  wave  amplitude  function  is  developed, 
and  three  methods  are  designed  for  the  estimation  of  a  ship’s  length.  The  first  section 
discusses  the  relationship  between  a  ship’s  hull  and  its  wave  amplitude  function;  the 
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second  section  gives  the  theoretical  proof  and  analysis  about  this  periodic  character 
for  general  ship  hull  shapes  and  the  relation  between  a  ship’s  length  and  its  bow  and 
stern’s  shapes;  the  final  section  gives  a  discussion  of  the  estimation  methods. 

5.1  Relationship  between  a  Ship’s  Hull  and  its  Wave  Am¬ 
plitude  Function 

In  the  study  of  fluid  motion,  wave  flows  due  to  a  moving  body  are  imagined  to 
be  generated  by  a  continuous  volume  distribution  of  singularities  within  the  body, 
extending  out  to  its  surface  [14].  In  this  section,  the  discussion  starts  directly  from 
the  relationship  between  the  wave  amplitude  function  and  wave  source  densities  or 
source  singularities.  The  non-dimensional  form  of  the  relation  has  been  described  by 
Eggers  et  al.  [20].  After  some  mathematical  manipulations,  the  dimensional  form  of 
this  relation  is  given  by 

AW  =  +  K„(0)V}dD  (5.1) 

where  a  represents  the  source  distribution,  D  is  the  source  region,  and  K ,  Kx,  and 
Ky  are  the  wave  numbers  defined  in  Chapter  2.  To  estimate  a  ship’s  length  and 
hull  by  using  this  formula,  the  following  two  assumptions  are  made  about  the  source 
region  and  the  relation  between  the  singularity  distribution  and  the  ship  hull. 

To  simplify  the  source  region,  it  is  now  assumed  that  the  ship  hull  is  thin,  that 
is,  the  beam  is  small  compared  to  all  other  characteristic  lengths  of  the  problem 
[14].  Thus,  the  singularity  distribution  can  be  envisaged  to  be  on  the  ship’s  center- 
plane,  instead  of  on  the  ship’s  hull  surface,  and  the  source  region  is  considered  to  be 
— %  <  x  <  %  and  —H  <  z  <  0,  where  L  and  H  denote  the  ship  length  and  draft, 
respectively.  Under  this  assumption,  (5.1)  becomes 

AW  =  —  J.  J_H°(x,z)exp{K{t)z  +  }K,(9)x)ixiz  . 


(5.2) 


This  formula  will  be  used  to  analyze  the  periodic  character  of  the  wave  amplitude 
function. 


In  the  estimation  of  ship  hull  shape  in  the  next  chapter,  the  explicit  relation 
between  the  wave  amplitude  function  and  the  ship  hull  is  useful,  and  thus  is  given 
below.  To  the  first-order  approximation,  the  relationship  between  a  hull’s  geometry 
and  its  centerplane  singularity  distribution  can  be  obtained: 

.  .  UdC(x,z) 

=  (5-3) 

where  £(x,  z)  defines  the  local  half-beam  of  the  hull  surface.  Thus,  combining  (5.2) 
and  (5.3)  gives  the  explicit  relation  between  the  wave  amplitude  function  and  the 
ship’s  hull: 

A (*)  =  ~  exP  {K(0)z  +  j  KT{0)x)dxdz  .  (5.4) 

For  simplicity  in  discussion,  the  normalization  of  x  and  z  with  respect  to  the  ship 
length  L  and  draft  H  is  considered.  By  letting 

X  =  |x'  x  €  [-1, 1]  (5.5) 

z  =  Hz'  z  €  [—1,0]  (5.6) 

equations  (5.2)  and  (5.4)  respectively  become 

A{0)  =  LHU3Kl(9 )  j\  fx  <r(x,  z')e^e)l'  e«9)*dxdz  (5.7) 

A{0)  =  ~  HUAKl{9)  J*x  J°x  9^\Z  emz'dxdz  (5.8) 


where 


V^°)  2 Kx^  2U*cos0  2F*  cos  0 

fi(0)  =  ¥H1k2(9)  =  — — —  =  — - - - 

^  )  g  (/2  cos2  0  L  F*  cos2  9 


(5.10) 


For  ease  of  notation,  the  prime  on  x  and  z  will  be  ignored  in  the  following  discussion. 
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5.2  Periodic  Character  of  the  Ship  Wave  Amplitude  Func¬ 
tion 

The  periodic  character  of  the  ship  wave  amplitude  function  can  be  observed  when 
the  wave  amplitude  function  is  described  as  a  function  of  the  longitudinal  wave  num¬ 
ber  Kx  and  plotted  on  a  Kx-A  diagram.  This  periodicity  is  proved  mathematically 
in  this  section,  and  the  inherent  connection  between  the  ship  length  and  the  periodic 
character  will  be  discussed. 

In  the  following  derivation,  the  wave  amplitude  function  with  two  dimensional 
integral  form  in  (5.7)  is  rewritten  into  a  marginal  integral  of  only  x.  Its  integrand 
is  expressed  in  the  form  of  a  power  series  and  the  integral  is  then  calculated.  Af¬ 
ter  mathematic  manipulations,  the  real  and  imaginary  parts  of  the  wave  amplitude 
function  are  expressed  in  the  form  of  cosine  functions  and  their  periodic  characters 
are  then  analyzed. 

The  wave  amplitude  function  in  (5.7)  can  be  expressed  as  a  marginal  integral  of 
x  by  defining  a  function  F(x),  i.e., 

A(Kk)  =  A{0)  =  J*  F(x)e^gdx  (5.11) 

where 

F(x)  =  ^LHU3Kl  1°  a(x,  z)  e^dz  .  (5.12) 

9  J- x 

Here,  the  wave  amplitude  function  has  also  been  written  into  a  function  of  variable 
Kx  instead  of  6.  The  integrand  function  F(x)  is  another  weighted  integral  of  the 
singularity  distributions.  With  the  assumption  given  in  (5.3),  F(x)  will  directly 
relate  to  the  ship  hull  shape,  thus  generally  F(x)  is  a  smooth  function.  In  the 
discussion  here,  it  is  assumed  that  F(x)  and  its  derivatives  are  continuous  in  the 
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region  —1  <  x  <  1,  so  that  F(x)  can  be  expressed  in  the  form  of  a  power  series  of  x 
F(x)  =  =  5Z  a2.*2‘  +  a2.+i*3<+1  •  (5.13) 

j=0  t=0  1=0 

Now,  substitute  (5.13)  into  (5.11)  and  then  integrate  it.  In  the  calculations,  the 
following  integral  formulas  Me  useful: 


Jxncosax  =  pRn(a,x)  cosax  +  ^r„(o,x)  sinax 
Jxn  sinax  =  pfn(a,x)  cos  ax  +  qin(a,x)  sinax 


(5.14) 

(5.15) 


where  n  is  an  integer  and 


PRn(a,x)  =  { 


for  n  =  0 


for  n  >  0 


[?1 


n! 


.n-2r 


9fln(<Z,x)  1)r(n-2r)!  a2^1 

m 


n! 


n-2r 


P/nKx)  E(  1)+1(n_2r)f  a2r+l 


?;,(«.*)  =  < 


0 


for  n  =  0 


(5.16) 

(5.17) 

(5.18) 

(5.19) 


l  ErJo  '(-li'irffcj)!  ■  £?Sf1  for  n  >  0 
With  the  above  integral  formulas,  the  wave  amplitude  function  in  (5.11)  becomes  a 
form  from  which  the  periodic  character  can  be  observed  much  more  easily: 


A(KX)  =  Ar(Kx)  +  )Aj{Kx) 

=  Qr{Kx)  cos  eR(Kx)  +  jQi(Kx)  cos  Q!(Kx) 


(5.20) 


where  Ar(Kx)  =  Qr(Kx)  cos  Qr(Kx)  and  Aj(Kx)  =  Q j{Kx)  cos  Qi{Kx )  denote  the 
real  and  imaginary  parts  of  the  wave  amplitude  function.  When  the  function  F(x ) 
is  an  even  function  of  x,  Ar  will  vanish.  With  the  notation  v  =  %KX  as  in  (5.9), 
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Qr(Kx),  Qi{Kx)  ,  Qr(Kz)  and  ©/(#*)  in  (5.20)  are  defined  by 


Qr{Kx)  =  2{(£a2,pflJi(i/))2  +  (£a2,gRsj(i/))2]*  (5.21) 

i=0  i=0 

Ql(K.)  =  2[  (£  aw  p„.„  (!.))=■  +  (f>„+1  („))’]*  (5.22) 

1=0  i=0 

6r(Kx)  =  V  -  <t>R{v)  (5.23) 

9i{Kx)  =  v-<i>i{v)  (5.24) 


where 


i  /..\  i-_—  lr  ^i=0  a2iQRii{,/)  i 

Mv)  -  ‘g&ai^W1 

«*)  =  tan~‘[  ] 

£,=0  ®2i+lP/a<+i  iy) 


prM  = 


for  i  =  0 


for  t  =  1,2,3,... 


=  S(~1)r(^_2r)» 

Ph« A”)  =  E(~1)r+l  (2/_  2r  +  1)!  ' 

qhi*^v)  =  (2*  -  2r)!  ‘ 


(5.25) 

(5.26) 


(5.27) 


(5.28) 


(5.29) 


(5.30) 


From  (5.20),  it  is  found  that  both  the  real  and  imaginary  parts  Ar  and  Ai  of 


the  wave  amplitude  function  consist  of  signals  with  an  /^-varying  magnitude  and 


an  FT*- varying  “frequency”.  That  is,  both  the  magnitudes,  Qr(Kx)  and  Q[(KX), 
and  the  instantaneous  “frequencies”,  ^  and  ±  change  with  Kx.  Here,  the 
frequency  concept  for  a  time  signal  is  used.  By  quoting  the  terminologies  in  telecom¬ 


munications,  the  real  and  imaginary  parts  look  like  two  signals,  both  amplitude  and 
angle  modulated  with  “carrier  frequency”  ■—  if  Kx  is  the  analog  of  time  t. 

In  the  above  discussion,  it  has  been  assumed  that  F(x)  and  its  derivatives  are 
smooth  and  continuous  in  the  hull  surface  region.  Thus,  there  is  one  dominant  fre- 
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quency  component  at  frequency  in  the  above  expressions  of  Ar  and  A/.  However, 
the  real  situation  may  not  be  so  perfect.  F(x)  may  be  a  piecewise  smooth  function, 
that  is,  there  are  some  discontinuous  function  and  derivative  points.  In  this  case, 
harmonic  frequency  components  will  appear.  For  example,  assume  that  there  is  one 
discontinuous  point  at  x  —  xb  6  (—1, 1),  and  F(x)  is  expressed  in  the  form 


F(x) 


|  Ejto  a3xj  f°r  -1  <  x  <  xb 


{  for  xfc  <  x  <  1 

Substituting  (5.31)  into  (5.11)  yields  the  results 


(5.31) 


Ar(Kx)  =  QR^K^cosii/ -  <j>Rl(v))  + 

Qr2(Kx,  xb)cos(xbv  -  <j>R2(v,xb))  (5.32) 

Ai(Kx )  =  Qh(Kx)cos(v  ~  <f>h(v))  + 

Qi2(Kz,  xb)  cos {xbu  -  <f>h (t>,  x6))  (5.33) 


where  the  magnitudes  Qr1,  Qr2,  Qix,  Qi2  and  the  phases  <£«,,  4>r2,  <t>ix,  <f>i2  are 
combined  with  the  coefficients  Pn,(u,  x),  qR,(v,  x),  p/t(u,x)  and  qi^v,  x),  defined  in 
(5.16)  —(5.19),  at  x  =  — l,xj,  or  1.  The  detailed  expressions  for  these  magnitudes 
and  phases  can  be  found  in  the  appendix.  In  (5.32)  and  (5.33),  cos(i/  —  ^r,(u))  and 
cos(i/  —  <f>ix(v))  represent  higher  frequency  components  generated  at  the  ship’s  bow 
and  stern,  i.e.,  x  =  ±1,  and  they  are  similar  to  those  in  (5.20)  with  frequency 
The  new  frequency  components  cos {xbv  —  <j>R2(v,xb))  and  cos {xbv  —  4>ix(v,xb))  are 
generated  by  the  discontinuous  point  at  x  =  xb,  and  their  frequency  is  lower  than 
—  since  |x&|  <  1.  When  xb  =  0,  Ar  and  A/  may  contain  a  direct  current  component. 

Generally  speaking,  each  discontinuous  point  in  the  function  F(x)  may  add  a 
new  frequency  component  to  Ar  and  A/,  and  the  frequency  is  always  lower  than 
~£.  As  mentioned  before,  the  function  F(x)  is  related  to  the  shape  of  a  ship  hull. 
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Thus,  the  increase  in  the  discontinuity  of  the  hull  surface  or  its  derivatives  will 
result  in  an  increase  in  the  discontinuity  of  F(x)  or  its  derivatives  and,  hence,  in 
an  increase  in  the  frequency  components  in  the  wave  amplitude  function.  According 
to  the  ship  wave  resistance  theory,  the  resistance  is  proportional  to  the  weighted 
integral  of  the  square  of  the  wave  amplitude  function  [14].  Hence,  this  increase  of 
the  frequency  components  may  result  in  an  increase  in  the  ship  wave  resistance. 
Therefore,  ship  hulls  are  usually  designed  to  be  smooth  so  that  the  resistance  can  be 
reduced.  Additionally,  note  that  the  coefficients  in  (5.16)  — (5.1&,  contain  the  factor 
xn-2r-i  Qr  xn-2r  jn  each  term,  and  are  small  for  x  —  xt  <  1  compared  to  those  for 
x  =  ±1.  Thus,  these  lower  frequency  components  generally  will  not  be  dominant 
as  found  in  real  examples.  Therefore,  the  following  discussion  will  still  focus  on  the 
problems  with  the  assumption  of  smooth  hulls. 

In  the  estimation  of  ship  length  from  Ar  or  Aj,  the  phases  Qr(Kx)  and  Qi(Kx) 
are  more  interesting  than  the  magnitudes  Qr(Kic)  and  Qi(Kk)  because  of  the  direct 
relation  of  the  ship  length  with  the  phases.  For  this  reason,  the  discussion  about  the 
estimation  of  ship  length  from  Ar  and  Aj  will  mainly  focus  on  phase.  It  is  found 
from  (5.23)  and  (5.24)  that  Qr{Kx)  =  | Kx  and  Qi{Kx)  =  ~KX  when  <f>R(v )  /) 
approaches  zero.  Thus,  the  curves  Ar{Kx)  and  Ai(Kx)  corresponding  to  the  real  or 
imaginary  parts  of  the  wave  amplitude  function  have  the  period  If  the  period  is 
measured,  then  the  ship  length  can  be  estimated. 

In  general,  however,  neither  <fo(i/)  nor  <f>i(v)  approaches  zero,  thus  the  estimation 
of  ship  length  becomes  complicated.  In  terms  of  (5.25)  and  (5.26),  the  phases  <f>R 
and  4>i  depend  not  only  on  v  but  also  on  the  coefficients  a,-  which  are  related  to  a 
ship’s  hull  shape.  The  relationship  between  the  phases  and  the  ship  shape  is  usually 
not  straightforward  and  obvious.  If  a  reasonable  approximation  is  made,  however, 
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more  insight  on  the  relation  can  be  still  gained. 

5.2.1  Approximation  of  A(KX)  for  larger  v 

If  a  large  value  of  v  is  assumed,  the  expressions  for  Qr,  Qi,  <f>R  and  d>i  can  be 
greatly  simplified  and  it  is  found  that  they  depend  on  the  function  values  and  the 
first  derivatives  of  F(x)  at  the  ship’s  bow  and  stern.  This  subsection  discusses  the 
validity  of  the  assumption  and  gives  the  approximate  expressions. 

With  the  assumption  of  large  u,  the  approximate  expressions  of  the  coefficients  in 
(5.27)  —(5.30)  are  given  first.  The  sums  of  pR2i,  qR2i,  9/,1+1  and  qi2i+1  in  (5.27)  —(5.30) 
consist  of  i  or  i  +  1  terms.  For  not  too  large  i,  these  sums  can  be  approximated  by 
their  first  or  second  term  when  v  is  large  enough.  With  this  approximation,  (5.27) 
—(5.30)  become 


prM 

.2  ... 

i-;  =  iPrM 

(5.34) 

\  oi") 

(5.35) 

Phi+Av) 

(5.36) 

as 

(2z  +  !)^  =  (2z  +  l)?/,  (i/) 

(5.37) 

for  i  =  0, 1,2, 3, ... 


The  errors  caused  by  the  above  approximation  are  dependent  on  variable  v  and 
index  i.  The  approximation  errors  for  i  <  7  are  plotted  in  Figures  5.1  —5.4  together 
with  the  curves  of  PR2(v),qRo(v),  P/,(v)  and  q/t( i/).  It  is  found  from  the  figures  that 
given  i,  there  is  a  value  such  that  the  approximation  errors  will  become  very  small 
when  v  is  larger  than  this  value.  For  instance,  the  errors  approach  zero  when  u  is 
much  larger  than  3  for  i  <  3  and  when  v  is  much  larger  than  7  for  i  <  7.  Recalling 
expression  (5.13),  the  largest  power  of  x  in  (5.13)  is  j  =  7,  corresponding  to  i  =  3, 
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and  j  =  15  to  i  =  7.  Since  a  ship’s  hull  shape  is  usually  smooth,  the  series  F(x) 
will  be  convergent  and  the  coefficient  a j  in  (5.13)  will  be  very  small  or  even  can  be 
neglected  for  large  j.  From  the  figures,  v  >  3  is  a  reasonable  assumption  for  most 
real  situations  where  i  <  3. 

In  addition  to  index  *,  the  range  of  v  needs  to  be  verified  to  make  sure  that  the 
approximation  errors  sue  small  enough.  To  understand  which  value  of  v  is  reasonable 
for  an  approximation,  v  is  here  related  to  the  Froude  number  since  the  scale  of  the 
Froude  number  is  available  in  ship  building  references.  By  the  definition  of  the 
Froude  number  given  in  Subsection  3.3.3,  it  follows  that 

u  = —K  -  L  9  1 

2  1  2  U2cos0  2F2cos0  ' 

Usually  0.1  <  Fn<  0.5  for  ships  [16];  thus,  v  ranges  from  to  ^5  correspondingly 

and  increases  with  the  increase  in  the  wave  angle  0.  From  (5.38),  the  minimum  v  is 

determined  by  umin  =  -—j.  Thus,  i/min  =  2  corresponds  to  Fn  =  0.50,  vmin  =  3  to 

Fn  =  0.41,  and  i /m,n  =  50  to  Fn  =  0.10.  Therefore,  the  assumption  that  v  >  3  is 

suitable  for  the  cases  where  F„  <  0.41,  in  particular,  for  most  merchant  ships,  for 

which  0.1  <Fn<  0.3  [16]. 

To  simplify  the  expression  of  the  phases  <f>R  and  <f>j,  substituting  the  approximate 
expressions  in  (5.34)— (5.37)  into  (5.25)  and  (5.26),  it  follows  that 


tan' 

—  tan-l[  —  ] 
v 


where  and  /?/  are  given  by 


0R 

01 


E”0  2iq* 

HZoa2i 

E”o(2»  +  l)q«+i 

a2i+l 


(5.39) 

(5.40) 


(5.41) 


(5.42) 
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By  using  the  triangle  relation  that  tan~x£  =  —  <an_1|,  where  “+”  is  taken  for 

positive  £  and  "  for  negative  £,  4>r{v)  can  be  expressed  in  the  form 

4>r{v)  =  ±-  +  <f>no(v)  (5.43) 

where 

4>rd{v)  =  -  tan_1[^]  .  (5.44) 

Now,  the  phases  6a  in  the  real  part  of  the  wave  amplitude  function  can  be  given 
from  (5.23)  and  (5.44)  by 

0a  (I<x)  =  i/-^Ao(i/)-(±|)  (5.45) 

In  terms  of  the  above  relations,  the  real  and  imaginary  parts  of  the  wave  amplitude 
function  can  be  now  rewritten  as 

Ar(Kx)  =  ±Qr(Kx)  sin(u  -4ro(v))  (5.46) 

MKX)  =  Qi(Kx)cos(v-Mv))  •  (5.47) 

The  signs  before  *  in  (5.45)  and  before  Qr  in  (5.46)  are  taken  as  before,  i.e.,  “+”  is 
taken  for  positive  @r  and  ”  for  negative  (3r. 

At  this  stage,  it  is  now  possible  to  relate  the  phase  parameters,  0r  and  /?/,  and 
the  magnitude,  Q(KX)  and  Qi(Kx),  to  the  function  F(x).  To  obtain  this  relation, 
break  up  F(x)  and  its  derivative  F'(x)  into  their  even  and  odd  parts,  i.e., 

F(x)  =  Fe(x)  +  F0(x)  =  a2ix2i  +  JT  a2i+1x2,+1  (5.48) 

i=0  i=0 

F\x)  =  F^x)  +  Fi(x)  =  f)  2iajiX2i~l  +  £(2 i  +  1  )ali+lx*  .  (5.49) 

i=0  i=0 

The  values  of  F(x)  and  F\x)  at  the  end  point  x  =  1  can  be  evaluated  now,  and  they 
are  given  by 

Fe(l)  =  f >2|  (5.50) 

i=0 
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W)  =  I>*+i  (5-51) 

isO 

K(  1)  =  f>a2i  (5.52) 

i=0 

^(1)  =  E(2.  +  I)«w  .  (5.53) 

1=0 

Comparing  the  above  expressions  with  the  definition  of  0r  and  0i  in  (5.41)  and 
(5.42)  leads  to 


0R 

01 


W) 
£(1) 
Foil)  * 


(5.54) 

(5.55) 


Note  that  the  even  and  odd  parts  of  a  function  can  be  expressed  in  the  function 
itself,  that  is, 


F.(x)  =  i[F(*) +  F(-*)1  (5.56) 

F.M  =  |[F(i)-F(-x)]  .  (5.57) 

Similarly,  the  derivative  of  F(x)  can  be  expressed  by 

F«'(x)  =  i(F'(x)-F'(-x))  (5.58) 

KM  =  i[F'(x)  +  F'(-x)]  .  (5.59) 


By  substituting  the  above  relations  into  (5.54)  and  (5.55),  0r  and  0i  can  now  be 
written  as 


F\  1)  -  F'(-l) 
F(l)  +  F(-l) 
F'(l)  +  F'(-1) 
F(l)-F(- 1) 


(5.60) 

(5.61) 


Similarly,  by  substituting  the  relations  in  (5.50)— (5.53)  and  (5.56)— (5.59)  into 
(5.21)  and  (5.22),  the  magnitudes  of  the  real  find  imaginary  parts  of  the  wave  am- 
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plitude  function  can  be  expressed  in  the  form 

Qr(K.)  =  i[(F'(l)-F'(-l))J  +  1/I(F(l)  +  F(-l))2]»  (5.62) 

Ql(K.)  =  i[(F,(l)  +  F'(-l)),  +  »*(F(l)-F(-l))*]i  .  (5.63) 

As  discussed  in  the  last  subsection,  the  variation  of  the  parameters  &r  and  /?/ 
causes  the  frequency  modulation  of  the  wave  amplitude  function,  and  complicates 
the  ship  length  estimation.  With  the  assumption  of  large  v ,  these  parameters  are 
simplified  and  are  related  to  only  the  values  of  the  function  F(x)  at  the  ship’s  bow 
and  stem.  Thus,  it  is  possible  to  find  their  values  by  solving  F(l),  F(— 1),  F'(l)  and 
F'(— 1)  from  (5.62)  and  (5.63).  Since  many  values  of  Qr  and  Qi  at  different  Kx  can 
be  obtained,  least  square  methods  can  be  used. 

So  far,  the  wave  amplitude  function’s  real  and  imaginary  parts,  Ar  and  A/,  have 
been  simplified  under  the  assumption  of  large  v .  From  (5.60)— (5.63),  it  is  found  that 
both  the  magnitude  and  the  phases  of  Ar  and  A/  depend  on  the  variable  v  =  \KX 
and  the  values  of  the  function  F(x)  at  the  ship’s  bow  and  stem.  This  result  also 
explains  the  phenomena  that  the  wave  generated  by  the  bow  and  stem  are  dominant 
in  the  ship’s  wake  compared  to  those  generated  by  other  parts  of  the  ship.  Thus, 
the  ship  wave  is  sometimes  considered  to  be  generated  by  a  moving  dipole,  or  a  pair 
of  moving  pressure  points  separated  by  a  distance  equal  to  the  ship  length.  Based 
on  the  approximate  expressions,  the  discussion  of  the  estimation  of  ship  length  from 
Ar  or  A/  is  further  given  in  the  following  subsections. 

5.2.2  Periodic  Character  and  the  Ship  Length  Estimation 

As  indicated  in  (5.46)  and  (5.47),  Ar  and  A/  represent  two  signals  that  are  both 
magnitude  and  angle  modulated.  Their  instantaneous  angle  frequencies,  and 
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are  given  from  (5.24)  and  (5.45)  as 


2m 


a  d  Qr  L 

“R  ~  ~dK~x  =  2  Z,*tf*  +  40& 

a  dQj  L  2L/3[ 

Uj  = 


(5.64) 

(5.65) 


dKx  2  L*K*  +  40}  ' 

These  two  expressions  show  that  the  angle  frequencies  have  a  direct  relation  with  the 
ship  length  and  wave  number  Kx.  In  order  to  understand  the  variation  of  the  angle 
frequencies  with  Kx  and  their  efFect  on  the  ship  length  estimation,  a  brief  analysis 
of  the  above  two  expressions  is  given  here.  Note  that  u>r  and  u;/  have  the  same  form 
except  for  the  subscripts.  Therefore,  the  discussion  below  will  focus  only  on  u>/,  and 
the  results  can  be  extended  to  u>r. 

For  the  extreme  case  where  Kx  approaches  infinity,  u ;/  approaches  a  constant  j-. 
By  determining  u>/  or  its  corresponding  period  7j,  the  ship  length  can  be  determined 
by 


Loo  =  2w/  =  ~  (5.66) 

where  L0 0  denotes  an  estimate  of  the  ship  length  when  Kx  approaches  infinity,  and 
it  depends  only  on  u/j  or  Tj. 

For  another  extreme  case  where  Kx  <C  u >/  approaches  another  constant 

~  fi)->  and  the  ship  length  can  be  determined  in  a  manner  similar  to  the  above 
case  if  w/  and  /?/  are  available.  If  0j  is  not  known  and  the  estimate  L «,  is  used  to 
determine  the  ship  length,  then  under-  or  over-estimation  may  occur  depending  on 
the  value  of  /?/,  i.e.,  positive  or  negative.  It  will  be  over-estimated  if  /?/  <  0  and 
under-estimated  if  /3j  >  0. 

For  general  cases,  the  wave  amplitude  function  is  recovered  from  wave  spectra; 
thus,  the  available  Kx  may  neither  approach  infinity  nor  satisfy  Kl  <C  In  these 
cases,  the  ship  length  is  determined  in  terms  of  (5.65).  If  is  used  to  determine 
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the  ship  length,  errors  will  be  introduced.  In  order  to  know  the  error  effect  on  the 
estimation,  consider  the  first  order  derivative  of  u>/.  From  (5.65), 

dui  4L3piKx 


U,  as 


(5.67) 


1  dKx  ( L*K l  +  4/3?  )2 

Because  Kx  >  0,  u'j  >  0  if  /3 /  >  0,  and  ui’j  <  0  if  /3/  <  0.  Hence,  the  angle  frequency 
u>i  increases  monotonically  with  Kx  when  fli  >  0  and  decreases  monotonically  with 
Kx  when  <  0.  That  is,  the  period  decreases  monotonically  with  Kx  when  /3/  >  0 
and  increases  monotonically  with  Kx  when  /?/  <  0.  Thus,  there  are  two  cases  for  the 
distribution  of  the  zero-crossing  points,  i.e.,  the  intersection  points  of  the  curve  A/ 
and  the  Kx- axis.  For  the  first  case,  where  /3/  >  0,  the  zero-crossing  points  becomes 


denser  as  Kx  increases;  for  the  second  case,  where  /3/  <  0,  the  zero-crossing  points 
becomes  less  dense  as  Kx  increases. 

According  to  the  above  analysis,  there  is  a  rule  of  thumb  to  know  the  trend  of 
the  estimation  error  when  is  used  to  determine  the  ship  length.  There  will  be 
an  over-estimation  if  the  zero-crossing  points  become  less  dense  (/3 j  <  0),  and  there 
will  be  a  under-estimation  if  the  cross-zero  points  become  denser  (/3 /  >  0).  This  rule 
is  also  suitable  for  the  estimation  from  the  curve  of  Ar{Kx).  The  above  estimation 
error  will  be  reduced  as  Kx  increases.  Hence,  it  is  suggested  that  the  period  Tj  at 
larger  Kx  be  taken  to  determine  L^.  In  Section  5.3,  several  ideas  will  be  proposed 
to  reduce  or  avoid  the  estimation  error. 

So  far,  the  periodic  character  and  its  effect  on  the  ship  length  estimation  have 
been  discussed.  The  conclusions  show  that  the  sign  of  the  value  of  fin  or  /3 /  plays 
an  important  role  in  the  estimation  performance,  and  the  sign  of  the  value  0r  or 
Pi  can  usually  be  determined  from  the  distribution  of  the  zero-crossing  points.  In 
some  cases,  however,  it  may  not  be  easy  to  identify  whether  the  zero-crossing  points 
become  more  or  less  dense.  For  these  cases,  the  phase  difference  between  the  curves 
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Ar(Kx)  and  Ai(Kx)  may  be  used  to  determine  the  sign  of  0r. 

To  determine  the  sign  of  (3r  from  the  phase  difference,  calculate  the  phase  differ¬ 
ence  from  (5.24)  and  (5.45): 

A0  =  0/  -  Qr  =  ±tj  +  (tan-1  -  tan-1  y^r)  .  (5.68) 

Note  here  that  before  y  “+”  is  taken  for  positive  0r  and  ”  for  negative  (3r.  In 
the  extreme  case  where  Kx  approaches  infinity,  the  value  of  (tan-1  jfj*-  —  tan-1  jfj*-) 
approaches  zero,  and  the  phase  difference  is  given  by  A0  =  *  if  Pr  >  0  and  A 0  =  — | 
if  Pr  <  0.  Thus,  the  sign  of  (3r  can  be  determined  according  to  whether  the  phase 
difference  is  -  or  — 

In  fact,  the  above  extreme  case  can  be  generalized  provided  that  the  value 
(tan-1  ~  tan-1  is  in  [  — y  ^  ].  If  the  phase  difference  is  limited  to  the  range 
(— ir,ir],  then  the  following  rule  of  thumb  to  determine  the  sign  of  /3r  is  obtained: 
0r  is  positive  if  0/  leads  Qr  (A0  >  0),  and  Pr  is  negative  if  0/  lags  behind  Qr 
(A 0  <  0). 

5.2.3  Periodic  Character  and  the  Shape  of  a  Ship’s  Bow  and  Stern 

This  subsection  reveals  a  relation  between  the  periodic  character  and  the  shape 
of  a  thin  ship’s  bow  and  stern  under  the  assumption  of  separation  of  variables.  This 
assumption  allows  the  ship  hull  surface  function  to  be  written  in  the  separated  form 

=  f(z)Hz)  .  (5.69) 

Note  that  here  x  and  z  are  the  normalized  variables.  The  function  £(x,  zo)  represents 
a  waterline  curve  of  a  ship  at  z  =  z0.  For  the  upper  half  waterline  curve,  C(x>^)  >s 
positive.  Thus,  positive  f(x)  and  h(z)  can  be  found.  When  C(xix)  is  separated  as 
in  (5.69),  the  waterline  curve  depends  mainly  on  the  function  f(x).  The  shape  of 
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the  curve,  particularly  at  the  bow  and  stern  of  a  ship,  has  an  effect  on  the  periodic 
character  of  the  wave  amplitude  function. 

With  the  above  separation  assumption,  the  phases  and  the  magnitude  of  An  and 
Ai  given  in  Subsection  5.2.1  can  be  simplified  and  related  to  the  hull’s  first  and 
second  order  derivatives.  For  a  thin  ship,  the  function  F(x )  defined  in  (5.12)  can  be 
written  in  the  following  form  by  substituting  (5.3)  and  (5.69)  into  (5.12)  : 

F(x)  =  f'(x)co(Kx)  (5.70) 


where 


co(Kx)  =  —HU^Kl  f°  h(z)e»*dz  . 
*9  J- i 


(5.71) 


Equation  (5.70)  is  substituted  into  (5.60)  —(5.63),  yielding  the  following  expres 
sions  for  the  phase  parameter  and  magnitudes. 

/"(I) -/"(-l) 


Pr  = 
Pi  = 


/'(l)  +  /'(-D 
/"(I)  +/"(-!) 


(5.72) 

(5.73) 


/'(l) -/'(-l) 

Qr(Kx)  =  i^[(/'(l)-r(-l))J  +  ^V(l)  +  /(-l))2l»  (5-74) 

Qi(Kr)  =  [(/"(!)  +  /"(-I))2  +  *2(/(l)-  /'(-l))2 ]*  (5-75) 

where  /* ( 1)  and  /*(—!)  denote  the  first  order  derivatives  at  the  bow  (i  =  1)  and 
stem  (z  =  —1),  and  /”(1)  and  /”(— 1)  denote  the  second  order  derivatives  at  the 
bow  and  stem. 

In  terms  of  the  function  f(x )  and  the  above  expressions,  the  relation  between  the 
parameters,  Pn  and  /?/,  and  the  geometric  shape  of  a  ship’s  bow  and  stern  can  be 
analyzed.  According  to  the  geometric  meaning  of  the  first  and  second  derivatives  of 
a  function,  the  slopes  /  (l)  and  /  (  — 1)  are  proportional  to  the  bow  and  stern’s  half 
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Figure  5.5:  Hull  waterline  curve  y  =  ±£(x,  z0). 

angles,  and  /(l)  <0  and  /  (  — 1)  >  0  for  the  upper  half  waterline  curve  as  shown  in 
Figure  5.5.  The  second  derivatives  /"( 1)  and  /"  ( — 1 )  describe  the  concaveness  and 
convexity  of  the  bow  and  stern.  When  the  second  derivative  at  the  bow  or  stern  is 
positive,  the  shape  of  the  waterline  curve  at  the  bow  or  stem  is  concave,  and  when 
the  second  derivative  at  bow  or  stern  is  negative,  the  shape  is  convex.  Therefore, 
the  parameters  and  /?/  depend  on  the  the  half  angle  and  the  concaveness  or 
convexity  of  the  bow  and  stern.  In  terms  of  (5.73),  if  both  bow  and  stern  are  convex 
(or  concave),  then  will  be  positive  (or  negative);  if  one  is  convex  and  the  other 
is  concave  the  sign  of  /?/  will  be  determined  by  whichever  shape  is  dominant,  i.e., 
whichever  shape  is  more  extreme.  A  similar  conclusion  can  not  easily  be  obtained 
for  from  (5.72),  because  the  values  of  /"(l)  —  /”(— 1)  and  f'(  1)  +  /(— 1)  depend 
on  the  specific  derivative  values. 
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The  above  analysis  tells  us  that  the  half  angles  and  the  concaveness  or  convexity  of 
the  bow  and  stern  are  related  to  the  parameters  0r  and  /?/.  Thus,  the  half  angles  and 
the  concaveness  or  convexity  of  the  bow  and  stem  can  be  predicted  if  the  parameters 
0R  and  are  known.  In  fact,  however,  the  sign  of  the  parameters  0r  and  /?/  cam 
also  be  used  to  predict  roughly  the  concaveness  or  convexity  of  the  bow  and  stem. 
Obtaining  the  sign  is  usually  much  easier  than  obtaining  the  exact  values  of  fa  and 
fa.  As  discussed  in  Subsection  5.2.2,  the  sign  of  fa  and  fa  may  be  found  from  the 
distribution  of  the  zero-crossing  points  and  the  phase  difference.  The  following  is 
the  main  conclusion  about  the  prediction  of  the  concaveness  or  convexity  of  the  bow 
and  stem  from  the  sign  of  fa  and  fa: 

1)  If  /?/  >  0,  then  the  bow  and/or  stern  have  a  convex  shape.  This  is  because 
/'(l)  +  /  (— 1)  <  0  when  fa  >  0;  thus,  f"(l)  and/or  /"(— 1)  are  negative. 

2)  If  /?/  <  0,  then  the  bow  and/or  stem  have  a  concave  shape.  This  is  because 
/  (I)  +  /”(—!)  >  0  when  fa  <  0;  thus,  /'( 1)  and/or  f(— 1)  (or  both)  are  positive. 

If  the  sign  of  the  parameter  fin  is  also  available  and  if  it  can  be  assumed  that  the 
half  angle  of  the  stem  is  larger  than  that  of  bow  or  vice  versa,  then  additional  infor¬ 
mation  about  the  concaveness  and  convexity  of  the  bow  and  stem  can  be  obtained 
by  analyzing  (5.72)  and  (5.73)  together. 

5.2.4  Ship  Length  Estimation  from  the  Magnitude  of  A{KX) 

In  the  above  subsections,  it  has  Seen  assumed  that  both  the  real  part  Ar(Kx)  and 
the  imaginary  part  A/(KX)  of  the  wave  amplitude  function  A(KX )  are  available.  In 
many  practical  situations,  however,  only  the  magnitude  of  the  wave  function  A(KX) 
is  available,  and  its  phase  is  not  known  or  ambiguous.  Thus,  Ar(Kx)  and  Ai(Kx )  can 
not  be  obtained.  This  particularly  happens  when  only  the  power  spectrum  is  known 
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and  the  phase  information  is  lost.  This  may  also  happen  in  the  case  where  the  ship 
center  can  not  be  exactly  determined.  This  is  because  when  the  wave  spectrum  is 
calculated  from  the  wave  elevation,  the  origin  of  the  wave  data  coordinate  system  has 
a  translation  with  the  ship  center,  i.e.,  the  origin  of  the  reference  coordinate  system; 
thus,  an  additional  phase  factor  is  produced.  Since  the  wave  amplitude  function  is 
recovered  from  the  wave  spectrum,  the  wave  amplitude  function  has  an  additional 
phase  factor,  too. 

This  subsection  discusses  the  estimation  of  ship  length  from  the  magnitude  \A(KX)\. 
The  periodic  character  can  also  be  found  in  |y4(A'x)|,  as  in  the  real  and  imaginary 
parts  of  A(KX),.  For  simplicity  in  theoretical  analysis,  the  square  of  the  magnitude, 
denoted  by  Am(Kz ),  is  considered  in  the  following.  From  (5.46)  and  (5.47),  it  follows 
that 

Am(Kx)  £  \A(KX)\2  =  po(Kx)  +  p(Kx)  cos(2i/  -  KW)  (5.76) 

where 

MK.)  =  QtiK^+QftK')  (5., -7) 

+<?)(*.)]*  (5.78) 

<4  (u)  =  tan- 1|  Q«(A'*)sin(2fe° Mi-Q/PQsinpM*')) 

1  '  l%(A-t)co»(2*R„(^))-<3i(A-I)cos(2*(^))  •  10 

Since  there  is  always  a  p  <  po,  the  right  hand  side  of  (5.76)  is  nonnegative  for  any 
Kx,  and,  thus,  the  absolute  symbol  is  omitted.  Am(Kx)  looks  also  like  a  signal,  both 
magnitude  and  angle  modulated,  but  the  ^arrier  frequency”  ^  is  double  compared 
to  that  of  Ar{Kx)  or  Ai(Kx),  The  frequency  again  contains  the  information 
about  the  ship  length  L. 
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To  analyze  the  periodic  character  of  Am(Kx),  the  magnitude  and  phase  axe  dis¬ 
cussed  further  here.  Comparing  Po(Kx)  with  p(Kx)co s(2v  —  4>m(v)),  it  follows  that 
po(Kx)  is  a  slowly  varying  component.  The  estimation  of  ship  length  is  based  on  the 
periodic  character  in  cos(2t/  —  <f>m(v)).  Hence,  we  need  to  make  sure  that  p(Kx)  is 
also  a  slowly  varying,  non-zero  component  compared  to  cos (LKX  —  4>m(v)).  To  do 
this,  the  phase  difference  (0/  —  0/*o)  in  (5.78)  )  is  first  simplified  below.  In  terms  of 
(5.40)  and  (5.44),  ( 0/  —  <f>no)  is  rewritten  as 

-  <f>R n{v)  =  tan~l[  -  ]«(£/-  Pr)~  (5.80) 

v  -  PiPr  v 

where  the  assumption  that  v2  >  02  and  v2  ft],  has  been  made  in  the  last  step. 
Now,  consider  the  behavior  of  p(Kx)  as  v  — ►  oo.  In  terms  of  (5.80),  cos  2(0/  —  <f>m) 
approach  1  when  v  becomes  very  large,  and  thus  the  oscillation  magnitude  of  p(Kx ) 
approaches  !((?//( if*)  —  Q2(KX))2 1.  By  using  the  results  in  Subsection  5.2.1,  p{Kx) 
can  be  expressed  for  large  v  from  (5.62)  and  (5.63)  as 


P(KX)  «  \Q2n(Kx)  -  Q2(Kx)\ 

=  (5.8i) 


Thus,  p(Kx)  changes  slowly  with  Kx  compared  to  cos (LKX  —  0m( v)).  If  v  >• 
I  ^(i)F(-~i)^  1»  t^ien  P(K*)  *s  approximately  equal  to  ^  \F(l)F(— 1)|. 

Now,  consider  the  phase  0m  for  large  v.  By  substituting  (5.40)  and  (5.44)  into 
(5.79),  0m  is  given  approximately  by 


<j>m(v)  «  tan~l[2v 


Q2M»2  +  0k)  ~  Q'rPr o(^2  +  0!) 


}  (5.82) 


Q2M2  ~  0f)(*  +  0k)  -  QU*  -  0k)(^  +  01) 

If  it  is  further  assumed  that  v 2  0)  and  v2  0\,  then  0m  can  be  approximated 

by 


4>m 


.  -iPm 
-tan  — 


v 


(5.83) 
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where 


QlPi  -  QlPm 
Q2i~Q2r 


Substituting  (5.60)  -(5.63)  in  (5.84)  yields 


lim  0„ 


F\l)  F\-l) 

F(l)  F(- 1)  * 


(5.84) 


(5.85) 


According  to  the  above  analysis,  the  phase  <j>m  will  be  small  when  v  is  large,  thus 
LKX  in  the  phase  plays  a  dominant  role  and  the  periodic  character  can  be  observed 
in  the  curve  Am(Kx).  Methods  for  estimating  ship  length  from  Am(Kx)  will  be 
introduced  in  Section  5.3.  One  practical  example,  estimating  the  length  of  a  real 
ship  hull  model  from  Am(Kx),  will  be  given  in  Section  5.4,  and  the  results  show  that 
good  estimation  can  be  obtained  from  Am(Kx). 


5.2.5  Examples  of  Ship  Hulls 

This  subsection  gives  some  examples  of  ship  hulls  to  demonstrate  their  wave 
amplitude  function  and  evaluate  the  parameter  /?  from  the  theoretical  calculation 
and  the  approximation  formulas  given  in  the  above  subsections. 

The  first  example  is  Wigley’s  Hull,  which  is  frequently  used  in  theoretical  analysis. 
The  normalized  expression  of  Wigley’s  hull  is  describe  by 

=  f  (i-z’K1-*2)  x  €  [-1,1],*  €  [-1,0]  (5.86) 

where  x  and  z  are  the  variables  normalized  by  the  half  ship  length  and  draft,  re¬ 
spectively,  as  defined  in  (5.5)  and  (5.6).  The  primes  have  been  omitted  for  ease  of 
notation.  Substituting  (5.86)  into  (5.8),  the  wave  amplitude  function  is  given  by 

A(K.)  = 


(5.87) 


where 


Qi(Kx)  =  ci(Kx)-—^y/A  +  L2K%  (5.88) 

Cl{Kx)  =  ThZ^kT)  [tl{Kx)2  ~ 2  + 2(1  +  ^Kx) )e~M(y,) ]  (5-89) 

H(KX)  =  —Kl  .  (5.90) 

9 

Indeed,  the  wave  amplitude  function  in  (5.87)  has  the  form  we  expect,  and  its  phase 
consists  two  terms.  The  second  term  in  the  phase,  — is  recognized  as  4>i 
as  discussed  before,  and  it  becomes  very  small  as  Kx  increase.  Thus,  the  first  term 
^  Kx  will  be  dominant. 

Now,  let  us  evaluate  the  parameter  j3j  from  the  theoretical  result  and  from  the 
approximation  method  directly  based  on  the  shape  of  hull.  Comparing  the  phase  in 
(5.87)  with  that  in  (5.24)  and  (5.40)  yields  =  1.  To  obtain  0i  directly  from  the 
hull  shape,  consider  the  hull  shape  function  C(*> z)-  Since  it  is  separable,  /?/  can  be 
directly  estimated  from  f(x)  —  1  —  x1  in  terms  of  (5.73),  and  has  the  same  value  as 
above.  Here,  /?/  is  larger  than  zero,  thus  the  zero-crossing  points  of  the  curve  Aj(Kx) 
become  denser  as  Kx  increases  and  under-estimation  may  occur. 

The  second  example  is  the  Cosine-Sine  Hull  with  the  normalized  hull  expression 

C(*,*)  =  +  cos(xx)][l  -sin(^z)]  x  €  (-1, 1],*  €  {-1,0]  (5.91) 

Similar  to  Wigley’s  hull,  its  wave  amplitude  function  can  be  found  and  is  given  by 

A(KX )  =  j  Qj(Kx)  sin(|  Kx)  (5.92) 

where 


Qi{Kx ) 


SirBU^Kxii{Kx) ,  1  f 

g(4x>  -  L'Kl)  1  ,<(**)  %(*,)’  + (f)2 


1  ,  »{Kt) 

»(KX)^  »{Kxy  +  ( f)’ 


)e^*>]  . 


(5.93) 
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This  theoretical  result  can  be  further  verified  from  the  analysis  of  the  ship  hull 
function.  For  this  hull,  f(x)  =  1  —  x3  and  /?/  approaches  infinity  when  it  is  evaluated 
in  terms  of  (5.73).  Thus,  «f>j  is  *  from  (5.40)  and  cos {%KX  —  <t>i )  =  sm(^Kx). 
This  is  the  same  as  the  above  theoretical  result  in  (5.93).  Note  that  the  angle 
frequency  is  a  constant  equal  to  y,  meaning  that  theoretically,  there  will  be  no  over- 
or  underestimation  of  ship  length  from  the  wave  amplitude  function. 

The  third  example  is  the  Wigley-Cosine  Hull.  Its  normalized  hull  expression  is 
described  by 

=  f(x)'{l~  z2)  x  6  [-1,1],*  €  [-1,0]  (5.94) 


where 


/(*) 


I  5  [1 +cos(;rx)]  0<x<l 


(5.95) 


[  1  —  x2  —  1  <  x  <  0 

This  hull  has  a  discontinuity  at  x  =  0,  i.e.,  the  derivatives  /"(0),/(4)(0), ...  are 
not  continuous.  Additionally,  it  is  not  symmetric  in  the  x-direction,  thus  the  wave 
amplitude  function  will  contain  both  real  and  imaginary  parts,  which  are  given  by 


where 


Ar{Kx)  =  Qm(Kx)  +  QR(KT)cos(-Kx-<j>R(v)) 

Ai(Kx )  =  QI{Kx)cos(^Kx-</>I{t')) 


(5.96) 

(5.97) 


Qro(Kx)  =  <*(*,)( 

Qr(Kx)  =  -Cl(tf.)(3  +  ( 
Qi{Kx) 

<MI/) 

<h{v) 


_ 1] 

2(^2  -  IT2)  ^  1 


tan-1{i/  [ 


+  —  )21$ 

v2  '  v2(»/2  —  it2)  T  v2’  J 

)2)* 


4,2  X2 


*  V  2(u2-ir2) 
4(i/2  —  x2) 


(4  +  ir2)i/2  —  Air2 


1} 


x2*/2 


—  tan_1{—  [1-  ■-  -  1} 

xu  1  4  {v2-ir2)u 


(5.98) 

(5.99) 

(5.100) 

(5.101) 


(5.102) 
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Here,  ci(Kx)  is  the  same  as  in  (5.89)  and  v  —  %KX.  To  compare  the  above  result  to 
that  from  the  approximation  method,  we  substitute  f(x)  =  1  —  x2  into  (5.72)  and 
(5.73),  and  obtain  the  parameters  0r  and  0i  in  (5.39)  and  (5.40),  i.e.,  0r  =  and 
0i  =  ,  which  «e  the  same  as  the  limits  of  the  factors  in  the  square  brackets 

in  (5.101)  and  (5.102)  as  v  — *  oo.  Since  0r  >  0,  the  estimation  from  Ar  will  be 
under-estimated,  and  since  0i  <  0  the  estimation  from  Aj  will  be  over-estimated  if 
no  compensation  is  made. 

5.3  Methods  of  Determination  of  Ship  Length 

As  discussed  in  previous  sections,  the  wave  amplitude  function  has  a  periodic 
character  and  the  ship  length  can  be  predicted  from  the  periodicity.  The  simplest 
method  to  predict  the  ship  length  is  to  evaluate  the  period  when  Kx  approaches 
infinity  and  then  calculate  as  described  in  Subsection  5.2.2.  Since  Kx  depends 
on  the  wave  angle  0 ,  i.e.,  Kx  =  a  large  Kx  requires  a  large  resolvable  wave 

angle.  In  real  situations,  however,  the  maximum  available  wave  angle  is  limited  by 
the  data  sampling  interval  and  ship  speed,  which  has  been  discussed  in  Section  3.3. 

In  general,  the  more  periods  available  in  the  data  of  A(KX),  the  better  for  ob¬ 
taining  a  good  estimation  of  ship  length.  However,  the  number  of  the  periods  or  the 
available  zero-crossing  points  depends  on  the  maximum  available  wave  angle.  The 
relation  between  the  wave  angle  and  the  number  of  zero-crossing  points  can  be  un¬ 
derstood  through  a  simplified  model  of  A(KX).  As  seen  in  the  above  discussion,  the 
real  or  imaginary  part  of  A(KX)  has  a  pattern  like  cos(^Kx  —  <f>),  the  zero-crossing 
points  appear  at 

~KX  —  <j>  =  ^(2m  +  1)  for  m  =  m0,m0  +  l,m0  +  2,  •  •  •  (5.103) 

and,  thus,  the  wave  angles  at  the  zero-crossing  points  can  be  expressed  from  (5.9) 
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and  (5.103)  in  the  form 


°  COS  1t((2m  +  l)x  +  *)*2I 


(5.104) 


where  Fn  =  is  the  Froude  number  and  mo  is  the  smallest  integer  such  that 
((2mTi)»4^)ii'1  —  Thus,  the  large  number  of  zero-crossing  points  needs  the  large 
maximum  available  wave  angles  available. 

Because  of  these  reasons,  Kx  can  not  be  very  large  in  practical  situations.  Hence, 
the  effect  of  the  phase  factors  4>r,  <j>j ,  or  <f>m  on  the  period  estimation  may  not  be 
negligible.  In  addition,  in  order  to  obtain  results  quickly  and  accurately,  data  pro¬ 
cessing  and  automatic  detection  are  necessary.  In  this  section,  several  algorithms 
to  determine  ship  length  are  introduced.  Although  the  discussion  about  these  algo¬ 
rithms  focuses  mainly  on  the  software  realization,  it  is  also  possible  to  use  them  in 
hardware  realizations  for  real  time  estimations. 


These  algorithms  include  the  spectrum  method,  zero-crossing  method,  and  fre¬ 
quency  demodulation  method.  In  the  spectrum  method,  the  period  is  estimated  by 
calculating  the  power  spectrum  of  Ar,  Aj  or  Am.  In  the  zero-crossing  method,  the 
zero-crossing  points  of  the  curve  Ar,  Aj  or  Am  are  detected,  and  then  used  to  find 
the  period  variation  with  Kx  for  further  ship  length  estimation.  In  the  frequency 
demodulation  method,  the  frequency  of  Ar,  Ai  or  Am  is  demodulated,  and  then 
the  frequency  variation  with  Kx  is  used  to  estimate  ship  length  .  In  the  following, 
these  methods  are  demonstrated  through  an  example  of  a  Wigley-Cosine  hull,  since 
its  wave  amplitude  function  contains  both  real  and  imaginary  parts  and  thus  is  a 
typical  example. 

Generally  speaking,  the  behavior  of  the  curves  Ar  and  Ai  is  better  than  Am,  thus 
the  estimation  of  ship  length  from  Ar  or  Ai  is  easier  than  from  Am.  However,  in 
practical  situations,  the  estimation  from  Am  may  be  more  useful  since  Am  is  easier 


Figure  5.6:  Spectrum  method  for  the  estimation  of  ship  length 
to  determine  than  Ar  or  A/. 

5.3.1  Spectrum  Method 

A  scheme  for  the  estimation  of  ship  length  using  the  spectrum  method  is  shown 
Figure  5.6.  In  general,  the  signal  Ar,  A/  or  Am  can  be  directly  inputted  to  the 
system  for  processing.  In  many  real  situations,  however,  the  available  data  record 
are  short,  only  a  couple  of  periods  of  signal  in  length,  in  particular,  for  Ar  and  A/. 
Thus,  it  is  preferred  to  use  Aj t  and  A]  as  the  input  signals,  so  that  the  frequency  of 
the  signals  is  doubled,  and  the  number  of  zero-crossing  points  increases.  This  will 
be  helpful  to  detect  the  signal  periods.  Before  the  square  operation  is  taken,  it  is 
suggested  that  the  signal  Ar  or  A/  be  filtered  using  a  high  pass  filter  to  remove  the 
direct  current  component  and  the  component  lower  than  as  shown  in  Figure  5.6. 
For  example,  consider  Ar(Kx)  =  Qrq(Kz)  +  Qr(Kx)cos(^Kx  -  4>r{v)).  After  the 
high  pass  filter,  the  lower  frequency  component  Qro(Kx)  is  filtered  out  and  the 
signal  becomes  Qr(Kx)cos(jKx  —  4>r(v)).  After  the  square  operation,  the  signal  is 
IQr(Kx)  [1  —  2 cos (LKX  +  2<^h)  },  where  cos (LKX  +  2 </>r)  is  the  desired  component 
and  has  frequency  If  the  high  pass  filter  is  not  used,  there  will  additionally  exist 
the  components  of  Qro{Kx)Qr{Kx)cos{^Kx  —  4>r)  which  is  not  desired  now.  After 
this,  the  squared  signal  will  undergo  further  processing  for  estimating  ship  length. 
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Hull  Waterline  Curves 


x(m) 

Figure  5.7:  A  Wigley-Cosine  hull  model  and  its  waterline  curves  y  =  db£(x,z)  at 
z  =  —0.3, —0.2, —0.1, 0.0  meters. 

In  the  bottom  digram  in  Figure  5.6,  a  bandpass  filter  centered  at  fam  —  £ 
is  used  to  remove  noise  and  undesired  lower  or  higher  frequency  components.  For 
example,  the  component  Qk(Ax)  in  the  squared  signal  is  filtered  out.  The  cutoff 
frequency  of  the  filter  can  be  roughly  determined  by  measuring  the  period  of  the 
input  signal  s,-.  Then,  the  output  signal  si(Kx)  of  the  filter  is  used  to  calculate 
the  power  spectrum.  At  this  stage,  the  spectrum  diagram  shows  a  pulse  at  around 
ffcz  —  If  the  transverse  axis  is  labeled  as  L  =  2 »■/*■,,  then  the  ship  length  can 
be  directly  read  from  the  pulse  position.  In  order  to  determine  the  pulse’s  frequency 
position  automatically  and  accurately,  however,  a  method  called  the  pulse  position 
detection  can  be  used.  In  this  method,  the  pulse  is  cut  by  a  threshold.  The  average 
frequency  position,  denoted  by  /* ^ ,  of  the  data  points  larger  than  the  threshold  is 
then  calculated  and  is  used  to  obtain  the  ship  length,  L^>  =  2x//cj0. 


Wave  Amplitidue 
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Figure  5.8:  (a)  Wave  amplitude  function  of  the  Wigley-Cosine  hull  model;  (b)  Fre¬ 
quency  variation  of  Ar  and  Aj  and  their  approximation. 


ooow  * 


Figure  5.9:  Signal  waveforms  in  the  ship  length  estimation  from  Ar  for  the  Wigley- 
Cosine  hull  model  using  the  spectrum  method.  The  estimated  ship  length 
is  L  —  4.825  m  with  a  relative  error  of  2.6%. 
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As  a  example,  consider  a  Wigley-Cosine  hull  model  with  length  L  =  4.953  meters, 
width  B  =  0.978  meters  and  draft  H  =  0.362  meters.  The  ship’s  speed  is  given  as 
U  =  2.229  meters /second.  The  hull  surface  and  the  waterline  contour  plot  are  given 
in  Figure  5.7.  The  real  part  Ar ,  the  imaginary  part  Aj  and  the  magnitude  |A|  of 
the  wave  amplitude  function  and  the  frequency  variation  of  Ar  and  A;  are  shown 
in  Figure  5.8.  Ar,  Ai  and  |A|  are  calculated  from  (5.96)  and  (5.97).  In  Figure  5.8, 


-  f  =  -ffli 

are  the  approximations  of  Ac or  and  Au>/,  and  are  calculated  from  (5.64)  and  (5.65) 


with  0r  =  and  0i  = 


The  signals  at  each  stage  in  Figure  5.6  are  gi\  <*n  in  Figure  5.9.  The  input  signal 
Si  is  the  squared  Ar,  the  real  part  of  the  wave  amplitude  function.  There  are  128 
data  points,  each  separated  by  an  interval  of  0.0445  rad./meter.  Kx  ranges  from 
1.974  to  7.626  rad./meter,  corresponding  to  the  wave  angle  from  0°  to  75°.  The 
cut-off  frequency  of  the  high  pass  filter  is  0.1  Hz,  and  the  cut-off  frequencies  of  the 
bandpass  filter  are  0.5  and  1.5  Hz.  The  spectrum  is  directly  calculated  according  to 
the  definition  given  in  [11]  using  the  FFT  algorithm.  The  spectrum,  s2,  has  been 
normalized  by  the  value  at  the  peak,  and  the  abscissa  has  labeled  directly  the  length 
scale  instead  of  2ir  Jkm  for  easy  observation  of  the  length  estimation.  Cyclic  spectral 
analysis  methods  and  other  advanced  spectral  analysis  methods  may  be  used  for 
high  performance  [35]  -  [38].  Here,  however,  a  simple  method  is  used  to  increase  the 
spectrum  resolution,  that  is,  the  data  length  is  increased  to  1024  points  by  padding 
zeros.  The  computer  calculation  gives  the  ship  length  estimation  a s  L  =  4.825 
meters.  It  is  under-estimated  as  predicted  from  theoretical  analysis,  with  a  relative 
error  of  2.6%. 
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Figure  5.10:  Zero-crossing  method  for  the  estimation  of  ship  length. 

5.3.2  Zero-Crossing  Method 

A  scheme  for  the  estimation  of  ship  length  using  the  zero-crossing  method  is 
shown  in  Figure  5.10.  Similar  to  the  spectrum  method,  the  input  signal  is  considered 
to  be  Am(Kx )  or  the  filtered  and  squared  signal  Ar{Kx)  or  A/(KX).  The  difference 
is  that  following  the  bandpass  filter,  the  signal  is  clipped  instead  of  being  used  to 
calculated  the  power  spectrum.  The  clipper  is  set  to  have  a  threshold  such  that  the 
signal  s2(Kx)  is  a  square  wave  signal  with  a  magnitude  close  to  zero.  Then,  the 
clipped  signal  is  differentiated,  giving  S3  =  The  differentiation  is  realized  by 
using  the  three  point  Lagrangian  interpolation  method  [34].  The  differentiated  signal 
is  combined  with  the  pulses  located  at  zero-crossing  points.  To  obtain  the  periods, 
the  positive  pulses  are  taken  using  another  clipper.  The  positive  pulse  signal  s4(Kx) 
is  then  used  to  calculate  the  periods  at  different  Kx.  Similarly,  the  pulse  position 
detection  method  is  used  here  for  fast  and  accurate  calculation.  The  interval  between 
two  pulses  is  the  period  at  the  corresponding  Ks.  The  first  and  last  pulses  may  not  be 
counted  to  avoid  false  periods  at  the  beginning  and  end  of  the  input  signal.  Finally, 
these  periods  are  used  to  calculated  the  ship  length. 

In  the  calculation  of  ship  length  from  the  periods,  different  strategies  can  be 
considered.  If  the  periods  {Txzi,i  =  1, ...,/}  do  not  form  a  monotonic  sequence 
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and  their  lengths  appears  to  be  random,  then  the  ship  length  is  calculated  from  the 
average  of  the  period  length,  that  is,  Lzc  =  YL\=i  Tk Mi-  This  case  may  happen 
in  the  estimation  from  Am(Kx),  or  the  magnitude  of  the  wave  amplitude  function. 
If  the  obtained  periods  {Txsi,i  =  1, ...,/}  do  form  a  monotonic  sequence,  then  the 
ship  length  is  calculated  from  the  last  period,  i.e.,  Lze  =  2irTfczl.  In  this  case, 
the  estimation  error  trends  can  be  predicted  in  terms  of  whether  the  sequence  is 
monotonically  increasing  or  decreasing.  If  the  variation  in  periods  is  small  near  TksI, 
then  the  error  will  not  be  large.  To  achieve  high  accuracy  in  the  estimation,  the 
following  method  may  be  used  if  at  least  three  zero-crossing  positions  are  available. 
For  the  signal  Am,  squared  Ar  or  A/,  the  phase  is  given  by  9  «  LKX  +  where 
/?  denotes  /?TO,  0r  or  /?/.  Since  the  cosine  function  has  a  phase  period  2n7r,  0  can  be 
written  in  the  form 

LK’i+~nh  =  i+2{i~l)*  i=1 . '  (5105) 

where  /  is  the  number  of  zero-crossing  points  to  be  considered.  If  three  or  more 
Kxi  are  available,  then  the  ship  length  and  (3  can  be  found  by  solving  this  set  of 
equations.  Least  square  methods  can  be  used  if  more  than  three  points  can  be 
available.  As  an  example,  a  three-point  formula  is  given  in  the  following.  Asstime 
that  there  are  three  points  Kxi,  Kx 2  and  Kx 3.  The  ship  length  and  f3  are  estimated 
from 

L  = 

fi  = 

The  above  discussion  considers  the  signal  s^AT*)  with  only  the  positive  pulses.  If 
the  negative  pulses  are  considered  too,  then  (5.105)  can  be  modified  as  follows. 

LKxi + ~nr  =  (2i "  x)\  *=1’~’r 


2ir 


(■ 


K- 


*3 


K.. 


xl 


Kx3  -  Kxl  xKx3-Kx2 
L2  Kx3  -  2Kx2  +  Kxl 

4  Arf  "  /fc  +  ^7 


Kx2  -  Kxi 
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Figure  5.12:  Frequency  demodulation  method  for  the  estimation  of  ship  length. 

where  /’  is  the  number  of  zero-crossing  points  Kx%  to  be  considered.  L  and  /?  can  be 
derived  in  the  same  way. 

As  an  example,  consider  the  same  Wigley- Cosine  hull  model.  Here,  the  input  is 
Aj.  The  signals  at  each  stage  in  Figure  5.10  are  given  in  Figure  5.11.  The  cut-off  fre¬ 
quencies  of  the  bandpass  filter  are  0.4  Hz  and  11.2  Hz  (Nyquist  frequency).  According 
to  the  calculated  results,  the  zero-crossing  point  position  is  3.55372, 4.79989, 6.04606, 
without  considering  the  first  and  last  points.  The  periods  are  1.24617  and  1.24617. 
Thus,  the  estimated  ship  length  is  L  =  5.042  meters.  It  is  over-estimated  as  predicted 
from  theoretical  analysis  with  a  relative  error  of  1.80%. 

5.3.3  Frequency  Demodulation  Method 

In  the  above  zero-crossing  method,  the  periods,  which  vary  with  Kx,  are  estimated 
by  locating  the  position  of  the  zero-crossing  points,  and  the  ship  length  is  then 
calculated.  In  the  frequency  demodulation  method,  the  instantaneous  frequency 
is  detected,  and  then  the  ship  length  is  estimated.  The  scheme  for  the  frequency 
demodulation  method  is  shown  in  Figure  5.12.  The  first  four  operations  in  the 
diagram  transform  the  input  signal  s,-,  an  unequal-amplitude  angle  modulated  signal, 
into  an  equal-amplitude  angle  modulated  signal;  the  last  four  operations  detect  the 
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“carrier  frequency”  jjj,  and  then  the  ship  length  is  calculated. 

In  order  to  obtain  an  equal-amplitude  signal,  a  clipper  is  used  following  the 
bandpass  filter.  Instead  of  determining  the  position  of  the  zero-crossing  points  as 
in  the  zero-crossing  method,  the  purpose  for  clipping  Si(Kx)  here  is  to  obtain  an 
equal-magnitude  square  signal.  Thus,  the  threshold  of  the  clipper  is  large,  compared 
to  that  in  the  zero-crossing  method.  Following  the  clipper,  a  bandpass  filter  is  then 
used  to  obtain  an  equal-amplitude  angle  modulated  signal.  After  the  amplitude  is 
normalized,  the  unit-amplitude  angle  demodulated  signal  S4  has  the  form 

s4(Kx)  =  cos(LKx  -  2 <f>)  (5.109) 

where  <f>  =  —  and  0  denotes  0m,  0r  or  0j.  The  “carrier  frequency”  ^  is 

desired  and  will  be  recovered. 

In  order  to  detect  the  frequency  information,  the  differentiation,  square  operation 
and  filtering  are  further  used.  The  output  of  the  differentiator  is  the  signal 


ss(Kx)  =  -(L--[r^-^)sin(LKx-2<f>)  .  (5.110) 

The  magnitude  of  contains  the  desired  frequency  information.  The  above  differ¬ 
entiation  can  be  realized  by  the  Fourier  transform  method.  The  method  is  based 
on  the  property  of  the  Fourier  transform.  That  is,  for  an  arbitrary  function  /(x), 
its  derivative  is  equal  to  the  inverse  Fourier  transform  of  jwjr{f(x)}.  A  discussion 
of  the  design  of  an  optimal  FIR  differentiator  can  be  found  in  [33].  To  obtain  the 
frequency  information  from  the  magnitude,  a  squarer  is  used  and  the  result  is 

MK.)  =  4  =  \(L  -  -  co<2LK’  -  W  ■  («•!”) 

Then,  a  low  pass  filter  is  used  to  filter  out  the  high  frequency  component  cos(2LKx  — 


K  (nd/m) 


7 


K  (radVm) 


Figure  5.13:  Signal  waveforms  in  the  ship  length  estimation  from  Am  for  the  Wigley- 
Cosine  hull  model  using  the  frequency  demodulation  method.  The  esti¬ 
mated  ship  length  is  Ljd  =  4.949  meters  with  a  relative  error  of  0.08%. 
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4<f>).  The  output  of  the  filter  is 

=  5(i- 15 ■  <5112) 

This  is  the  desired  frequency  information  containing  the  ship  length.  Similar  to  the 
zero-crossing  method,  several  strategies  can  be  used  to  estimate  the  ship  length.  If 
s7(Kx)  is  not  a  monotonic  function,  then  0  is  not  a  constant  and  the  average  value 
of  s7(Kx)  can  be  considered.  Since  usually  L2K2  4|/3  —  02\,  the  estimate  of  ship 
length  is  given  as 

L  »  y/2f7  (5.113) 

where  s7  is  the  average  of  the  signal  s7(Kx).  If  s7(Kx)  is  a  monotonic  or  almost 
monotonic  function,  then  0  can  be  considered  as  a  constant  and  the  ship  length  can 
be  obtained  by  solving  (5.112).  Since  there  are  two  unknowns,  L  and  0,  in  (5.112) 
and  many  function  values  of  s7(Kx)  are  available,  least  square  methods  can  be  used 
here. 

As  an  example,  consider  the  Wigley-Cosine  hull  model  again.  Here,  the  input  is 
Am.  The  signals  at  each  stage  in  Figure  5.12  are  given  in  Figure  5.13.  The  cut-off 
frequencies  of  the  filters  are  0.4  Hz  and  1.5  Hz  for  the  first  bandpass  filter,  and  0.5 
Hz  and  1.5  Hz  for  the  second  bandpass  filter.  The  cut-off  frequency  for  the  low  pass 
filter  is  0.03  Hz.  The  signal  s7  in  the  diagram  represents  the  frequency  of  Am,  the 
ordinate  has  been  directly  labeled  with  the  length  scale  L  =  2t/kx  instead  of  /ks 
for  easy  observation  of  the  estimated  length.  The  dotted  line  is  the  average  of  the 
frequency,  which  gives  the  estimation  of  ship  length.  The  estimated  ship  length  is 

A 

Lfm  =  4.949  meters  with  a  relative  error  of  0.08%.  The  standard  deviation  of  the 
variation  of  2ir/*,  is  0.166  meters. 


100 


So  far,  three  methods  for  the  estimation  of  ship  length  have  been  introduced.  In 
general,  the  spectrum  method  is  a  fundamental  and  powerful  method,  in  particular 
for  signals  with  noise  or  other  undesired  frequency  components.  The  zero-crossing 
method  and  frequency  demodulation  method  may  be  used  to  achieve  more  accurate 
results.  In  real  situations,  the  methods  can  be  combined  in  use.  For  example,  the 
spectrum  method  may  be  used  first  to  obtain  a  rough  estimation,  then  the  result  can 
be  refined  using  the  zero-crossing  method  or  the  frequency  demodulation  method. 

5.4  Estimation  of  the  Quapaw  Hull  Length 

This  section  applies  the  above  theoretical  analysis  and  methodology  to  a  practical 
problem,  the  estimation  of  the  Quapaw’s  length.  This  ship  model  has  length  L  = 
4.953  meters,  width  B  =  0.978  meters  and  draft  H  =  0.362  meters.  Here,  the 
assumptions  that  B/f  <1  and  B/H  is  small  are  not  valid;  thus,  this  is  not  a  thin 
ship,  that  is,  the  hull  is  not  narrow  and  deep.  The  hull  surface  and  the  waterline 
contour  plot  axe  given  in  Figure  5.14.  The  hull  was  towed  in  a  tow  tank  with  a 
constant  speed  U  =  2.229  meters/second  (7.31  feet/second).  Three  wave  elevation 
cuts  were  measured  with  three  independent  sensors  fixed  in  the  tow  tank  as  the  ship 
hull  passed.  These  wave  cuts,  shown  in  Figure  4.9  in  Chapter  4,  are  parallel  to  the 
ship  centerline  and  at  the  distance  of  y\  =  1.219,  yj  =  1.524  and  y$  =  1.8288  meters 
from  the  centerline,  respectively.  The  wave  data  have  been  preprocessed  to  remove 
the  pulse-type  noise  caused  by  instrumentation  and  the  stationary  waves  caused  by 
the  finite  length  of  the  tow  tank. 

Generally,  four  steps  can  be  taken  for  estimating  the  ship  length,  that  is,  1) 
calculate  the  FFT  of  the  wave  elevation  rj(x,yi)’,  2)  calculate  the  wave  amplitude 
function  A(KX);  3)  roughly  guess  the  ship  length  from  the  curve  A(KX)\  4)  estimate 
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Hull  Waterline  Curves 


x(m) 

Figure  5.14:  Quapaw  hull  model  and  its  waterline  curves  y  —  :t£(x, z)  at  z  —  —0.3, 
—0.2,  —0.1,  —0.001  meters. 

the  ship  length  using  the  methods  introduced  above.  Figure  5.15  gives  the  FFT  of 
|-4(/^r)|  and  the  spectrum  of  detrended  |A(A'r)|  for  the  data  RUN3. 
Analysis  of  the  plots  in  Figure  5.15  is  given  now.  In  the  figure,  (a)  plots  the 
magnitude  of  the  FFT  of  the  wave  amplitude  function;  (b)  plots  the  part  of  (a) 
to  see  the  detail  for  lower  spatial  frequency  components,  and  the  abscissa  has  been 
relabeled  as  the  wave  number  Kx,  which  is  equal  to  2 xu  here.  It  is  found  from 
(b)  that  there  are  two  peaks  which  are  much  larger  than  others.  The  first  peak 
in  Figure  5.15(b),  at  Kx  =  1.974  rad. /meter,  corresponds  to  the  transverse  waves; 
the  second  peak,  at  Kx  «  2.5  rad. /meter,  corresponds  to  the  diverging  waves.  This 
second  peak  and  is  larger,  and,  thus,  the  wave  component  it  represents  is  dominant  in 
the  wave  cut.  The  dominant  wavelength  can  be  estimated  to  be  A  =  ~  =  2.5  meters. 
According  to  the  theoretical  analysis  in  Chapter  3,  the  wave  number  in  the  x-direction 


102 


u  (1/m)  Kx(i*d7m) 


Figure  5.15:  (a)  Magnitude  of  the  FFT  of  the  wave  elevation  for  data  RUN3-C; 

(b)  magnification  of  plot  (a);  (c)  the  magnitude  of  the  wave  amplitude 
function  |./4(if*)|  calculated  from  (a);  (d)  the  spectrum  of  detrended 

\A(Km)  |. 
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is  Kx  —  The  transverse  wave  crests  are  almost  perpendicular  to  the 

ship  centerline,  i.e.,  9  «  0;  the  transverse  wave  number  is  Kq  =  1.974  rad./meter. 
Since  the  wave  cut  is  parallel  to  the  the  ship  centerline,  the  apparent  transverse 
wavelength  in  the  wave  cut  is  the  same  as  Kq.  However,  the  apparent  diverging 
wavelength  in  the  wave  cut  is  larger  than  the  diverging  wavelength  itself  because  the 
wave  cut  is  not  perpendicular  to  the  diverging  wave  crests. 

An  empirical  formula  to  guess  the  apparent  diverging  wavelength  in  wave  cuts  or 
to  guess  ship  length  may  be  obtained  from  a  simplified  model  of  the  wave  amplitude 
function.  To  do  this,  the  relation  between  the  wave  amplitude  function  and  the  FFT 
of  the  wave  elevation  cut  given  in  (3.50)  is  rewritten  by  changing  variable  9  into  Kx 
in  the  form 

A(KX )  =  H&-)  —  JU'Kl  -  g2  (5.114) 

L“K  irg  v 

where  #(f*)  is  the  Fourier  transform  of  the  wave  cut  »?(x,y0).  When  U*K%  »  g2, 
the  magnitude  of  A(KX )  is  approximated  as 

|A(Jf,)i  »  Iff (§01  (5.115) 

or  the  magnitude  of  is  written  as 

|ff(|5)|  «  .  (5.116) 

Now,  consider  a  hull  symmetrical  in  length.  Its  wave  amplitude  function  can  be 
described  as  a  very  simple  form,  A(KX)  =  j Qcos[%(Kx  -  Kq)  +  4>k0],  where  <f>K0  is 
the  phase  at  Kx  =  Kq  and  Q  is  a  constant.  A(KX)  approaches  zero  when  Kx  <  Kq 
since  the  minimum  wave  number  is  K0.  Although  A(KX)  has  a  cosine  form,  \H(^)\ 
decays  very  fast  as  Kx  increases  because  of  the  factor  Thus,  the  first  maximum 
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as  Ks  >  Kq  is  important,  and  it  is  approximately  given  as 


K,m.,  «  K0  + 


(5.117) 


where  k  can  be  between  0.5  and  1.0,  depending  on  the  phase  <t>x0.  For  example,  k 
can  be  taken  as  1.0  for  </>k0  «  (2n  +  1)|  and  0.5  for  </>k0  ~  (2n  +  1)*.  From  this 
empirical  formula,  the  apparent  wavelength  in  a  wave  cut  can  be  roughly  estimated, 
if  the  ship  length  is  known,  using 


A„  = 


2x 


2ir 


K 


Kq  +  Kj 


2  *f;l 


(5.118) 

(5.119) 


1  +  kvF* 

where  Fn  is  the  Froude  number.  The  ship  length  can  be  roughly  estimated  using 


L  ~  K 


**_  -  Kq 


(5.120) 


From  the  FFT  of  the  wave  elevation  cut  in  Figure  5.15,  for  example,  KXmtx  «  2.5, 
thus  the  ship  length  is  roughly  estimated  as  3.14  to  6.28  meters. 

Figure  5.16  (c)  is  the  magnitude  \A(KX)\  of  the  wave  amplitude  function.  It  is 
not  so  smooth  as  the  one  in  the  above  example  of  the  Wigley-Cosine  hull  model. 
However,  the  peaks  labeled  with  P\  and  P2  still  can  be  recognized  as  a  period,  which 
is  approximately  equal  to  1.5  rad./meter.  Thus,  the  ship  length  can  be  roughly 
estimated  as  L  «  =  4.2  meters. 

To  analyze  the  frequency  components  of  |A(A*)|  in  detail,  the  data  |v4(jFiTx)|  is 
detrended  and  then  its  power  spectrum  is  calculated.  In  order  to  be  consistent  with 
the  ship  length  scale,  the  abscissa  has  been  labeled  as  2 k/k,  instead  of  the 
frequency  corresponding  to  Kx,  in  the  spectrum  diagrams  given  below.  According 
to  the  theoretical  analysis  in  the  above  sections,  Ar(Kx)  or  Ai(Kx)  contains  the 
component  cos(%Kx  —  (j>),  and,  thus,  has  a  dominant  peak  at  in  its  spectrum; 
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\A(KX)\  contains  the  component  cos(LKx—<f>),  and,  thus,  has  a  dominant  peak  at  L  in 
its  spectrum.  All  the  analysis  is  based  on  the  linearized  free-surface  condition  and  the 
thin  ship  assumption.  In  the  spectra  of  |A(ATX)|  calculated  from  real  ship  wave  data, 
however,  there  may  exist  many  peaks  presenting  different  frequency  components, 
as  seen  in  Figure  5.15  (d).  In  this  diagram,  the  peak  at  around  5  meters  is  the 
second  harmonic  component  L,  representing  the  ship  length.  The  basic  component 
at  j  is  a  small  peak  here  since  this  spectrum  is  of  |y4(/C*)j.  In  addition  to  these  two 
components,  there  are  several  higher  order  harmonic  components  and  other  frequency 
components.  These  additional  components  may  result  from  the  nonlinear  effects  of 
ship  waves,  which  have  been  neglected  in  the  above  theoretical  analysis.  These 
nonlinear  effects  on  the  wave  elevation  are  very  small,  and  wave  components  caused 
by  the  nonlinear  effects  have  a  very  small  energy.  However,  they  locate  at  the  end  of 
higher  frequencies  in  the  wave  elevation  spectrum;  thus,  these  components  are  greatly 
enhanced  in  \A(KX)\  because  of  the  factor  K x  in  \A{KX)\  «  This  is 

the  situation  found  in  Figure  5.15  (d),  in  which  a  couple  of  frequency  components 
are  larger  than  the  component  at  2 x/k,  =  5.  Thus,  we  need  to  guess  the  ship 
length  first  using  the  first  two  important  peaks  in  the  diagram  of  \A{KX)\  and  using 
the  empirical  formula  (5.120)  to  know  which  component  is  our  desired  one  in  the 
spectrum  of  |A(ATX)|. 

The  results  of  ship  length  estimated  form  the  three  wave  cuts  of  the  data  RUN3 
are  shown  in  Figures  5.16,  5.17,  and  5.18.  The  magnitude  |-A(/fx)|  of  the  wave 
amplitude  function  is  first  calculated  from  the  FFT  of  the  wave  cut,  and  then  it 
is  detrended.  The  detrended  signal,  denoted  as  Ad(Kx)>  is  filtered  by  a  bandpass 
filter  with  bandwidth  0.4  Hz  and  centred  frequency  0.8  Hz.  Finally,  the  filtered 
signal,  denoted  as  Af(Kx),  is  used  to  estimate  the  ship  length.  The  three  methods 


methods 

wave  cut  A 

wave  cut  B 

wave  cut  C 

3-cut  average 

rms  error 

Ltp  (cl) 

4.746(4.2%) 

4.812(2.8%) 

5.347(8.0%) 

4.968(0.3%) 

0.330(6.7%) 

i>xc  (^l) 

5.240(5.8%) 

5.190(4.8%) 

5.290(6.8%) 

5.240(5.8%) 

0.291(5.9%) 

Lfd  (cl) 

4.959(0.1%) 

4.946(0.1%) 

5.078(2.5%) 

4.994(0.8%) 

0.084(1.7%) 

Table  5.1:  Quapaw  hull  length  estimated  from  the  three  wave  cuts  of  RUN3  using 
the  spectrum  method,  zero-crossing  method  and  frequency  demodulation 
method.  The  unit  of  the  ship  length  in  the  table  is  meters ,  and  the  true 
hull  length  is  4.953  meters. 

introduced  in  the  last  section  are  used  here.  The  spectrum  of  A/  has  been  normalized 
in  the  figures.  The  results  of  the  ship  length  estimation  is  listed  in  Table  5.1,  where 
Lap,  Lzc  and  L/d  denote  the  estimated  lengths  using  the  spectrum  method,  zero¬ 
crossing  method  and  frequency  demodulation  method,  respectively,  and  ti  denotes 
their  relative  error.  According  t'  the  statistical  analysis  of  the  nine  lengths  estimated 
from  the  three  wave  cuts  using  the  three  methods,  the  mean  of  the  estimated  ship 
length  is  5.068  meters  with  a  relative  error  of  2.3%.  The  root  mean  square  (rms)  is 
0.243  meters  with  a  relative  rms  error  of  4.9%.  This  example  shows  that  although 
the  theoretical  analysis  and  the  methodology  are  based  on  the  linearized  free-surface 
condition  and  the  thin-ship  assumption,  they  can  be  used  for  estimating  a  non-thin 
ship’s  length  from  real  wave  data,  and  good  estimation  results  can  be  obtained. 


Figure  5.16:  Quapaw  hull  length  estimation  from  the  tow  tank  wave  cut  RUN3-A 
at  y  =  1.219  m  using  three  methods.  The  estimated  ship  length  is 
=  4.746  m  (  cr,  =  4.2%  )  for  the  spectrum  method,  LIC  =  5.240  m 
(ti  =  5.8%)  for  the  zero-crossing  method,  and  L/d  =  4.959  m  (ti  = 
0.1%)  for  the  frequency  demodulation  method. 


wave  elevation  (cm) 


Figure  5.17:  Quapaw  hull  length  estimation  from  the  tow  tank  wave  cut  RUN3-B 

at  y  =  1.524  m  using  three  methods.  The  estimated  ship  length  is  I 

Lrp  —  4.812  m  (cl  =  2.8%)  for  the  spectrum  method,  Lxc  —  5.190  m  ■ 

(ti  =  4.8%)  for  the  zero-crossing  method,  and  Lfd  =  4.946  m  = 

0.1%)  for  the  frequency  demodulation  method.  ■ 
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Figure  5.18:  Quapaw  hull  length  estimation  from  the  tow  tank  wave  cut  RUN3-C 
at  y  =  1.8288  m  using  three  methods.  The  estimated  ship  length  is 
L,p  =  5.347  m  (ej,  =  8.0%)  for  the  spectrum  method,  Lzc  =  5.290  m 
(ec  =  6.8%)  for  the  zero-crossing  method,  and  L/j  =  5.078  m  {cl  = 
2.5%)  for  the  frequency  demodulation  method. 


CHAPTER  VI 

EXTRACTION  OF  SHIP  HULL  GEOMETRY 

INFORMATION 


The  prediction  of  a  ship’s  hull  geometry  from  the  ship  generated  wave  pattern 
can  be  considered  as  an  inverse  Kelvin  wake  problem.  Solutions  to  this  problem  are 
becoming  more  practical  with  advances  in  remote  sensing  technology. 

The  idea  used  here  to  extract  a  ship’s  hull  shape  is  based  on  the  relation  of  the 
ship  hull  shape  and  the  wave  amplitude  function  under  the  thin-ship  assumption, 
and  the  relation  of  the  wave  amplitude  function  and  the  ship  wave  spectrum.  These 
relations  have  been  discussed  in  the  previous  chapters,  thus,  this  chapter  focuses 
on  the  technique  to  extract  a  ship’s  hull  geometry  shape  from  the  wave  amplitude 
function.  Mathematically,  the  problem  is  to  find  the  function  of  hull  surface  £(x,z), 
given  the  wave  amplitude  function  A(9)  and  the  integral 

A' (6)  =  H  f1  1°  dC(?’*)ejW9)*eM(g)«^z  (6.1) 

J-i J-i  ox 

where  v  and  n  have  been  given  in  (5.9)  and  (5.10),  and  A!{9)  is  defined  as 

jjl 

A'(0)  =  cos3(0)A(0)  .  (6.2) 

2 9 

In  the  above,  (6.1)  is  obtained  from  (5.8)  under  the  thin-ship  assumption,  x  and  z 
are  the  variables  normalized  by  the  ship  length  and  draft,  and  eM^*  is  a  kernel 

function  with  parameter  9. 
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To  solve  this  continuous  inversion  problem,  we  prefer  to  convert  it  to  a  discrete 
inversion  problem  first,  and  then  solve  it  using  discrete  inversion  techniques.  The 
reasons  for  this  are  that  it  is  difficult,  in  general,  to  solve  a  continuous  inversion 
problem  analytically  and  that  the  values  of  the  wave  amplitude  function  are  usually 
available  in  discrete  form.  To  obtain  the  discrete  form,  the  integral  may  be  approxi¬ 
mated  as  a  summation  using  the  trapezoidal  rule  or  some  other  quadrature  formula. 
In  this  hull  inversion  problem,  however,  the  spectral  method  is  adopted  instead  of 
the  above  quadrature  methods.  In  the  spectrad  method,  the  unknown  function  is 
approximated  as  a  weighted  summation  of  some  basis  functions.  Because  the  basis 
functions  aure  known  and  the  integrals  associated  with  them  can  be  computed,  solv¬ 
ing  the  continuous  inversion  problem  becomes  a  matter  of  solving  a  set  of  discrete 
equations  to  determine  the  weights  or  coefficients.  The  advantage  of  the  spectral 
method  is  that  the  number  of  unknowns  to  be  determined  can  be  reduced  greatly. 
This  is  especially  important  for  two  dimensional  inverse  problems. 

When  both  the  magnitude  and  phase  of  the  wave  amplitude  function  are  well 
known,  the  ship  hull  inverse  problem  becomes  a  lineair  inverse  problem.  The  singu¬ 
larity  of  the  kernel  matrix  in  the  linear  problem  and  noise  of  the  observed  data  may 
have  severe  effects  on  the  inversion  results,  and  thus  speciail  techniques  may  have 
to  be  considered.  In  particular,  the  methods  of  constrained  linear  inversion  and  the 
maximum  likelihood  estimation  with  constraints  will  be  considered  in  solving  the 
ship  hull  inversion. 

However,  if  only  the  magnitude  of  the  wave  amplitude  function  is  available,  then 
a  complicated  non-linear  inverse  problem  must  be  handled.  For  the  hull  inverse 
problem,  the  non-linear  inversion  is  solved  based  on  the  criterion  of  the  maximum 
likelihood  estimation  with  constraints  using  optimization  techniques.  Examples  of 
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linear  inversion  and  non-linear  inversion  will  be  demonstrated  for  both  mathemati¬ 
cally  well-defined  ship  hulls  and  the  model  of  the  real  ship  Quapaw. 

In  the  following,  bold  lowercase  letters  denote  vectors,  and  bold  uppercase  letters 
denote  matrices. 

6.1  Spectral  method 

In  this  section,  the  spectral  method  is  used  to  convert  the  continuous  inverse 
problem  to  a  discrete  inverse  problem.  In  the  following,  the  solution  of  the  inverse 
problem,  £(x,  z),  is  approximated  by  a  weighted  summation  of  basis  functions.  Insert¬ 
ing  the  summation  into  the  integral  equation  (6.1)  yields  a  set  of  algebraic  equations 
for  different  values  of  the  parameter  0.  The  solution  of  the  algebraic  equations  gives 
the  weights,  or  coefficients,  of  the  summation,  and  any  value  of  £(x,  2)  in  the  specified 
region  can  be  evaluated  from  the  summation  of  basis  functions. 

For  this  goal,  let  the  basis  functions  in  the  x—  and  z— directions  be  &(x),i  = 
1,2, and  =  1,2 respectively.  C(x,z)  IS  now  expressed  in  the 

form 

c(*^)  =  5m°o^(;r)^(2)  •  (6-3) 

<'=i  7=1 

and  its  partial  derivative  with  respect  to  x  is  given  as 

=  EE  «,*.(*)/*;(*)  (6.4) 

ax  i=l j= 1 

where  4>i{x)  denotes  the  derivative  of  <£,(x).  Substitution  of  (6.4)  into  (6.1)  gives  the 
following  linear  equation  with  the  coefficients  to  be  determined: 

A'(»)  =  E E  .«[  WKi{«)  +  iWtiiW ]  , 

i=l  ;=i 


(6.5) 


where 


WRii{9)  =  WXRi{0)WZj{0)  (6.6) 

WUj(Q)  =  WXIi(9)WZj(0)  (6.7) 

=  J  ^  Mx)  cos  [v(d)x)dx  (6.8) 

Wjr/i(0)  =  J  Ux)mnWW*  (6.9) 

WZj(9)  =  H  J°Jj(z)e^dz  (6.10) 


In  general,  the  amplitude  function  A(9)  is  complex,  and,  thus,  A'(0)  is  too.  If 
A'r(9 )  =  i?e{A'(0)}  and  A!t(0)  =  Im{A'(9)}  denote  the  real  and  imaginary  parts, 
respectively,  then  (6.5)  gives 

M  N 

4w  =  EE*i't'»»(s) 

'=i i=i 

A',W  ='£'£<niW„iW  .  (6.ii) 

t=l  )Sj 

Given  K  values  of  9,  there  is  a  set  of  equations 

iml  jml 

4(0*)  =  EX> 

i*l  ;=1 

lk  =  .  (6.12) 


In  order  to  get  a  unique  solution  of  the  M  x  N  coefficients,  a^  ,  there  must  be  at  least 
M  x  N  equations.  Thus,  the  number  of  parameters  0k.  should  satisfy  K  >\M  x  N . 

The  above  equations  can  be  written  in  a  vector  and  matrix  form  by  rearranging 
the  subscripts.  According  to  the  rule  that  the  subscripts  (({tj},j  =  1, 2,  •  •  • ,  N),  i  = 
1,2, are  mapped  to  ({/},/  =  \,2,-  •  •  ,MN),  and  the  coefficients  {(avj,j  = 
1,2, ,N), i  =  1,2, are  mapped  into  {6/,/  =  1,2, •••  ,MN},  (6.12)  can  be 
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written  in  the  form 

a  =  Wb  (6.13) 

where  a  =  [A*(0i),  A'r(92),  •  •  ■ ,  A'r(0k ),  A'j(0 i),  A!j{0 2),  •  •• ,  A'^Ok)]7  with  dimension 
2 K  x  1,  b  =  bMii]T  with  dimension  MN  x  1,  and  T  signifies  transpose. 

The  matrix  W  may  be  called  the  kernel  matrix  as  it  is  related  to  the  kernel  function. 
It  has  dimension  2 K  x  MN  and  is  given  as 


Wr m 

Wjuji 

Wtuifi 

•••  Wj?(v-i)Ari 

Wrmu 

•  •  •  WRMN\  ■ 

"W 

Wrim 

Wrm 

•••  V-1)N2 

Wruk 

WRtjH2 

W/Ul/C 

Wrmk 

WRlffK 

Wjai/r 

•••  Wr(m-\)NK 

WrmiK 

■  ■  ■  WRMNK 

W'mx  • 

■  M'/ijvi 

wnn 

Winn 

■  ■  ■  WIMS1 

JV/113 

Wn  w 

Wjim 

W'iju 

••• 

Wit4\1 

■  •  •  WIMN2 

■  W/llA- 

hw  • 

WI1NK 

WW 

...  Wj(U-l)NK 

M'/Vilf 

•  •  •  WIhiNK  - 

where  Wmjk  and  W7«i*  denote  Wfuj(0k)  Wuj{0k ),  respectively. 

Now,  the  continuous  inverse  problem  has  been  converted  to  an  algebraic  inverse 
problem.  In  (6.13),  a  is  usually  obtained  from  measurements,  W  is  known  from 
construction,  and  b  is  to  be  solved  for.  The  kernel  matrix  W  depends  on  the  wave 
angle,  ship’s  length,  draft  and  speed,  and  the  choice  of  the  basis  functions.  Thus, 
the  next  step  is  to  select  a  set  of  appropriate  basis  functions. 


6.2  Selection  of  Basis  Functions 


The  most  preferable  basis  functions  should  have  the  properties  of  easy  compu¬ 
tation,  completeness  and  rapid  convergence  [40].  The  property  of  easy  computa¬ 
tion  means  not  only  that  the  basis  functions  themselves  should  be  easily  calcu¬ 
lated,  but  also  that  they  should  make  the  integrals  in  the  elements  of  the  above 
matrix, WxRi{Q),WxRi{8)  and  Wzj(0),  be  analytically  integrable.  Thus,  costly  nu¬ 
merical  computation  time  can  be  avoided.  The  completeness  of  the  basis  functions 
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make  it  possible  for  the  solution  to  be  represented  to  an  arbitrarily  high  degree  of 
accuracy  by  increasing  the  number  of  basis  functions.  The  property  of  rapid  con¬ 
vergence  allows  the  number  of  basis  functions  used  to  be  as  small  as  possible.  In 
addition,  the  geometry  of  the  problem  is  another  fact  that  needs  to  be  considered 
when  the  basis  functions  are  selected.  Study  on  spectral  methods  indicates  that 
the  best  choice  for  95%  of  all  applications  is  an  ordinary  Fourier  series  or  a  set  of 
Chebyshev  polynomials.  For  the  non-periodic  problem,  the  Chebyshev  polynomials 
are  extremely  robust  and  give  good  results  in  almost  all  situations  [40]. 

Based  on  the  above  principles,  the  Chebyshev  polynomials  Tn(0  are  chosen  in 
both  the  x—  and  z— directions.  The  explicit  expressions  of  the  Chebyshev  polyno¬ 
mials  jT„(£)  are  given  as 

To(Z)  =  1 
Ti{$  =  t 

m)  =  2?-\ 

•  •  • 

WO  =  W-uo  (6.15) 

*€[-1,1];  n>  1. 

They  have  two  main  properties:  1)  Tn(0  is  even  for  even  n,  and  odd  for  odd  n;  2) 
Tn(il)  =  1  for  even  n,  and  T„(±l)  =  ±1  for  odd  n. 

The  basis  functions  <f>i(x)  and  /3j(z)  are  constituted  by  the  Chebyshev  polynomi¬ 
als.  One  strategy  to  select  the  basis  functions  is  to  set 


M*)  =  Ti(x )  i  *'  =  1,2 ,-  -,M 
fc(»)  =  W*)  ,  3  =  1.2, •••,1V. 


(6.16) 

(6.17) 
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The  reason  not  to  include  the  lowest  order  of  the  Chebyshev  polynomials,  T0(x),  in 
the  basis  functions  of  the  x— direction  is  that  the  derivative  of  TQ(x)  =  1  is  equal  to 
zero,  and  thus  whether  To(x)  is  included  or  not  has  no  effect  on  in  the  integral 
(6.1).  Even  if  Tq(x)  is  included,  its  coefficient  can  not  be  determined  by  solving 
(6.13).  Thus,  a  constant  must  be  finally  determined  according  to  the  boundary 
condition  of  £(x,  z)  after  the  coefficients  a,-;  are  found. 

In  fact,  the  function  representing  the  ship  hull  shape,  £(x,z),  is  zero  on  the  ship 
hull  boundary  or  outside  of  the  boundary  on  the  x  —  z  plane,  that  is,  in  the  normalized 
coordinate  system, 


<(±1,2)  =  0,  (6.18) 

<(x,-l)  =  0.  (6.19) 

These  conditions  can  be  imposed  on  the  basis  functions  with  the  advantage  that  no 
constant  needs  to  be  determined.  According  to  this  idea,  the  basis  functions  in  the 
x— direction  are  set  to  be 

f  7<+1(x)  -  1  if  *  is  odd  ; 

<f>i(x)  =  4 

I  7i+i(x)  —  x  if  i  is  even  , 

i  =  (6.20) 


If  it  is  desired  that  the  boundary  condition  in  the  z— direction  given  in  (6.17)  be 
imposed  on  the  basis  functions,  then  they  should  be  set  to  be 

Tj(z)  +  1  if  j  is  odd  ; 

Tj(z)  —  1  if  j  is  even  , 

3  —  1,2,  •  •  • ,  Af.  (6.21) 


AW  = 


The  proper  choice  of  basis  functions  may  not  only  make  the  elements  of  the  kernel 
matrix  W  easier  to  compute,  but  also  make  the  discrete  inverse  problem  simpler.  In 
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solving  the  inverse  problem  given  in  (6.1),  the  basis  functions  will  be  chosen  from 
(6.17)  and  (6.20).  The  basis  functions  given  in  (6.20)  have  the  property  that  <t>i(x) 
is  even  for  odd  i,  and  odd  for  even  i  according  to  the  even  and  odd  property  of  the 
Chebyshev  polynomials.  Thus,  <j>i(x)  is  even  for  even  i,  and  odd  for  odd  i.  This  fact 
further  results  in  that  Wxm{6)  =  0  for  odd  t  and  Wxn(6)  =  0  for  even  t  since  the 
integrands  in  the  symmetric  integrals  (6.8)  and  (6.9)  are  odd.  These  properties  of 
Wjuj{6)  and  Wnj(0)  can  be  used  to  simplify  the  kernel  matrix  W  in  (6.14).  To  show 
this,  consider  an  example  with  M  =  4,  N  =  2  and  K  =  4.  In  this  case,  the  matrix 
W  is  given  as 


0  0  VV/QJ!  Wf02 1  o  0  B'jMii  W/J421 

0  0  VVjrju  Wjtm  0  0  *W  Wrm  a 

0  0  Wja  13  Wfoaa  0  0  Wr43  3 

0  0  Wj}214  WjU2  4  0  0  WRtlt  Wr 424 

W/l3i  o  o  Wj3  u  W/jj,  o  o 

VV/m  W/1M  0  0  VV/jij  W/saa  0  0 

WhM  ^IlM  0  0  W,313  W/323  0  0 

.  W/U4  W/1J4  0  0  W/3,4  W/324  0  0 


(6.22) 


Observing  the  pattern  of  the  above  matrix,  we  can  find  that  if  the  columns  of  the 


matrix  are  rearranged,  then  the  matrix  can  be  written  in  the  form 


WjBU  WR321  Wjmii  Wjun 

w«13  Wr222  WR*12  WRA22 

WjB13  Wjq 33  ^Vr413  Wfl423 

WfUU  WfOH  W/ttU  Wr 4J4 

0  0  0  0 

0  0  0  0 

0  0  0  0 

.  0  0  0  0 


0  0  0  0 

0  0  0  0 

0  0  0  0 

0000 
Wnn  wii2i  ^/sn  W/sai 

W'/na  WVim  ^312  W'rsaa 

W'/llS  W'/m  ^313  ^333 

^JIH  W/124  W/314  W/3J4  . 


(6.23) 


The  new  pattern  of  the  matrix  shows  that  the  linear  equation  given  in  (6.13)  can 


be  split  into  two  independent  sets  of  linear  equations  with  smaller  dimensions  if  the 


matrix  elements  are  rearranged  according  to  the  rules  that  for  Wmjk,  the  subscripts 
{({*i}»j  =  1>2, •••,7V),»  =  2, 4,---, M'}  are  mapped  into  ({/i},/i  =  1,2, •  •  •  ,LX), 
where  M'  =  M  if  M  is  even,  and  M‘  —  M  —  1  if  M  odd,  and  where  L\  =  ^*iV; 
that  for  Wiijk,  the  subscripts  {{{ij},j  —  1, 2,  •  •  • ,  N),i  =  1, 3,  ■  •  • ,  M"}  are  mapped 
into  ({kJjfj  =  1,2,  •  •  •  ,£2),  where  M"  =  M  —  1  if  M  is  even,  and  M"  =  M  if  M 

n 

odd,  and  where  Z,2  =  ^3- IV.  In  addition,  the  coefficients  {(o,j,i  =  1,2, = 
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2,4,- •  •  ,M'}  are  represented  by  {cj,/  =  1, 2,  •  •  • ,  Li),  and  the  coefficients  {(a,j,j  = 
1, 2,  •  •  • ,  N),  i  =  1, 3,  •  •  • ,  M")  are  represented  by  {d/,  /  =  1,2,-*-,  Li).  After  the  sub¬ 
scripts  are  rearranged,  the  two  independent  sets  of  linear  equations  can  be  expressed 
in  the  matrix  forms 


where 


a*  =  W*c  (6.24) 

a,  =  W/d  (6.25) 

a„  =  (6.26) 

a,  =  [Ai(91),A',(6J),---,4(«*))T  (6.27) 

c  =  [c„cJ,...,c£,l]r  (6.28) 

d  =  [d1,d7,---,dL,]T  .  (6.29) 


The  matrix  W«  with  dimension  K  x  L\  and  W /  with  dimension  K  x  are  given 
by 


r 

wKai 

•  •  • 

WR*ll 

Wm*i'-2)N  i 

wiu#'n 

^rw'ni 

W*  = 

Wiai* 

...  W/Bwa 

WRM 

ia 

WRM'  N2 

(6.30) 

.  wmiK 

**r« »K 

•  •  •  WR3NK 

WMXK 

wrm‘  NK  . 

and 


w,= 


f  wnn 
W'/m 


Wnn 

wim 


wnNi  Wr3U  ■  W/(v"-a)m  wim"u 

VV/3U  •••  Wj  ^im"  la 


W 


m"  m 


JM"  W3 


« 6.31) 


L  WH1K  WntK  •••  W/SlX  •••  Wr(kt"-2)NK  WIU"lK  W1U"SK  j 


Using  the  above  mapping  rule  of  the  coefficients  a,y,  the  hull  surface  function  can 


be  expressed  in  a  concise  form 


C(x,*)  =bZ(x,z)c  +  h*{x,z)d 


(6.32) 
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where  vectors  he(x,z)  and  h0(x,  z)  are  defined  as 
hj(x,z)  = 

Mx)Pi(z),  •••  tM'(x)Pi(z)’  •••  4>m'(x)Pn(z)}  (6.33) 

h'(x,z)  =  [&(x)A(z),  Mx)Mz)i  -Mx)Pn{z), 

Mx)0l(z)i  -<i>M"(x)Pi(z)i  •'■<f>Mtt(x)pN(z)\  •  (6.34) 

In  (6.32),  hf(x,  z) c  represents  the  odd  part  and  bj(x,  z)d  the  even  part  with  respect 
to  x  if  the  basis  functions  $,(x)  are  selected  to  be  even  for  odd  i  and  odd  for  even  i 
as  in  (6.20). 

So  far,  two  sets  of  linear  equations  have  been  established  for  determining  the 
vectors  c  and  d,  and  thus  the  coefficients  a,j.  If  a  ship  hull  is  assumed  to  be  bow-to- 
stern  symmetric  in  its  longitudinal  distribution  of  volume,  then  £ (x,z)  =  C{~xi  z), 
and  thus  c  =  0,  and  =  0  since  the  real  part  of  the  wave  amplitude  function  is 
equal  to  zero.  In  this  special  case,  only  (6.25)  needs  to  be  solved. 

For  the  basis  functions  composed  of  the  Chebyshev  polynomials,  both  <j>i{x)  and 
Pi(z)  are  polynomials  of  x  or  z,  thus  Wxm,  Wxn  and  WZj  in  (6.8)-(6.10)  are  com¬ 
binations  of  integrals  like 


==  J  xn  cos[i/(0)x)]  dx 

(6.35) 

=  J  xn  sin[t/(0)x)]  dx 

(6.36) 

=  J°^  zne^zdz  . 

(6.37) 

All  of  these  are  explicitly  integrable,  and  Cn  =  0  for  odd  n  and  Sn  =  0  for  even  n. 
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6.3  Effects  of  Singularity  and  Data  Error  on  Linear  Inver¬ 
sion 

Theoretically,  once  the  measurement  data  a/j  and  a/  are  given,  and  the  elements 
of  the  matrices  Wr  and  W /  are  calculated  for  each  given  value  of  parameter  9 ,  the 
unknown  vectors  c  and  d  can  be  found  by  solving  (6.24)  and  (6.25).  To  obtain  a 
correct  or  reasonably  accurate  solution,  however,  we  may  have  to  take  account  of  the 
ill-condition  of  the  matrices,  the  noise  involved  in  measurement  data  or  the  both.  In 
the  situations  of  ill-conditioned  matrices,  the  solution  may  be  unstable  when  there 
is  a  very  small  noise,  or  even  no  noise,  just  for  a  computer’s  finite  floating  point 
precision.  Therefore,  some  special  measures  must  be  taken  to  solve  linear  inverse 
problems. 

6.3.1  Effects  of  the  Ill-condition  of  Matrices 

The  singularity  or  ill-condition  of  a  matrix  can  be  formally  defined.  According 
to  the  linear  algebra  theorem  [42],  any  M  x  N  matrix  A,  whose  number  of  rows  M  is 
greater  than  or  equal  to  its  number  of  columns  N,  can  be  written  as  the  product  of 
an  M  x  N  column-orthogonal  matrix  U,  an  N  x  N  diagonal  matrix  A  =  [diag(A,-)]  with 
non-negative  elements,  and  the  transpose  of  an  N  x  N  orthogonal  matrix  V,  that  is, 

A  =  U  (diag(Ai)]  VT  .  (6.38) 

If  the  matrix  A  is  square,  N  x  N  say,  then  the  inverse  of  A  is  given  by 

A-1  =  V  [diag(Y~)]  UT  .  (6.39) 

Ai 

The  condition  number  of  a  matrix  is  defined  as  the  ratio  of  the  largest  of  the  Aj’s 
to  the  smallest  of  the  Aj’s.  A  matrix  is  singular  if  its  condition  number  is  infinite, 
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Figure  6.1:  Elements  A,-  of  the  diagonal  matrices  decomposed  from  (a)  W r  and  (b) 
W /  with  K  =  135  and  Lx  =  6  associated  with  M  =  4  and  N  =  3;  (c) 
W/i  and  (d)  W /  with  K  =  135  and  L\  =  18  associated  with  M  =  6  and 
6. 

and  it  is  ill-conditioned  if  its  condition  number  is  too  large,  that  is,  if  its  reciprocal 
approaches  the  computer  machine’s  floating  point  precision. 

Using  this  definition  to  examine  the  matrices  W*  and  Wj  in  the  above  ship  hull 
inverse  problem  described  by  (6.24)  and  (6.25),  we  unfortunately  find  that  even  if 
there  are  no  errors,  the  matrices  may  become  ill-conditioned  even  for  not  very  large 
M  and  N.  As  an  example,  consider  the  matrices  W/j  and  W /  with  the  basis  functions 
given  in  (6.17)  and  (6.20)  and  with  parameters  L  =  100  meters,  H  =  10  meters  rind 
U  =  10  m/s  (equivalently,  the  Froude  number  F„  =  0.32  ).  Figure  6.1  illustrates  the 
elements  A,-  of  the  diagonal  matrices  decomposed  from  Wr  and  W /  with  parameters 
K  =  135  and  Lx  =  6  associated  with  M  =  4  and  N  =  3,  and  with  parameters 
K  =  135  and  L\  —  18  associated  with  M  =  6  and  N  =  6,  respectively.  From  the 
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Figure  6.2:  Kernel  function  e*W*  with  Fn  =  0.32  and  ^  =  0.10  . 

figure,  it  can  be  seen  that  the  last  elements  of  the  diagonal  matrices  are  close  to 
zero  and  their  reciprocal  approaches  a  very  large  number,  especially  for  the  matrices 
with  larger  dimensions.  The  condition  numbers  are  calculated  as  7.642  x  104  and 
6.333  x  104  corresponding  to  and  W /,  respectively,  for  M  =  4  and  N  =  3,  and 
1.018  x  1016  and  3.084  x  1017  corresponding  to  Wr  and  W /,  respectively,  for  M  =  6 
and  N  =  6.  Thus,  the  matrices  in  the  latter  case  are  considered  to  be  ill-conditioned. 
The  reason  that  the  matrices  become  ill-conditioned  for  the  large  number  of  basis 
functions  is  that  the  number  of  rows  K  in  the  matrices  W r  and  W /  increase  as  M 
and  N  increase  since  there  must  be  K  >  ^-N  for  even  M  or  K  >  for  odd  M 

to  avoid  an  underdetermined  solution,  but  one  row’s  elements  axe  not  very  different 
from  another  row’s  as  the  number  of  rows  increases.  For  example,  let  us  examine 
the  elements  of  matrix  Wj  in  (6.31).  From  (6.7),  (6.9)  and  (6.10),  the  elements  can 
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be  written  as 

WIijk  £  Wwfa)  =  (-  £  *,(x)  sin  \v{9k)x\dx)  •  (tf  £  ^(zjc^cte)  (6.40) 

Note  that  the  kernel  in  the  second  integral  varies  little  with  the  wave  angle 

0 ,  especially  when  9  is  in  0°  —  40°,  as  shown  in  Figure  6.2.  It  is  the  smoothness  of 
the  kernel  that  makes  the  rows  of  the  matrix  (6.31)  become  almost  linear  dependent. 
From  the  viewpoint  of  matrix  theory,  the  rank  of  the  matrix  is  decreased.  Therefore, 
an  increase  in  the  number  basis  functions  does  not  improve,  but  may  worsen  the  ill- 
condition  of  the  inverse  problem.  To  avoid  an  ill-conditioned  problem,  the  number  of 
basis  functions  should  be  reasonable,  although  a  good  approximation  for  the  function 
of  the  solution  may  be  achieved  for  a  large  number  of  basis  functions. 

The  condition  numbers  Nc  of  the  kernel  matrices  W*  and  W /  depend  on  the 
ratio  H/L  and  the  Froude  number  Fn  since  their  elements  depend  on  v  and  fi,  which 
are  related  to  H/L  and  the  Froude  number  F„,  as  seen  above.  Figure  6.3  shows 
the  relations  between  Ne  and  Hf  L  and  between  Nc  and  Fn  for  the  given  numbers  of 
basis  functions  M  =  4  and  N  =  3.  The  values  of  9  are  evenly-spaced  in  the  range 
[5°,80°],  and  K  —  135.  From  the  figure,  the  curves  of  Nc  versus  H/L  are  flat;  Nc's 
have  smaller  values  when  H/L  is  larger  than  0.08  or  when  Fn  ranges  from  0.2  to  0.5. 
These  areas  with  lower  condition  numbers  are  typical  in  real  situations. 

6.3.2  Effects  of  Data  Noise 

In  practical  problems,  the  amplitude  function  A{6)  =  AR{9)-\-jAn(9)  is  obtained 
from  direct  or  indirect  measurements;  thus  there  must  be  some  systematic  or  random 
noise  or  errors  associated  with  them  and  these  errors  further  cause  errors  in  the 
solution  c  and  d  of  (6.24)  and  (6.25).  The  solution  c  and  d  will  be  considered 
unstable  if  the  errors  in  them  are  unacceptably  large.  In  addition,  note  that  all  the 
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elements  of  the  matrices  Wr  and  W /  are  calculated  from  the  integrals  (6.8)-(6.10) 
which  contain  the  parameters  of  the  ship’s  length  L  and  draft  H.  The  errors  of  these 
two  parameters  will  cause  errors  in  the  matrices,  and  finally  increase  the  errors  in 
the  solution.  This  subsection  will  focus  on  the  error  analysis  of  the  solution  c  and  d. 

First,  measured  or  observed  data  are  expressed  into  their  true  values  plus  noise. 
It  is  assumed  that  the  noise  on  the  measured  data  Ajt(0k)  and  >1/(0*)  is  additive 
noise,  thus,  >1/1(0*)  and  >1/(0*)  with  noise,  denoted  as  j4r(0*)  and  >1/(0*),  can  be 
expressed  as  sums  of  the  true  values,  Ar(0*)  and  >1/(0*),  and  noise  AAr(0*)  and 
AA/(0*): 

Ar(0*)  =  >1r(0*)  4-  A>1r(0*) 

i/(0fc)  =  A/(0*)  +  A>i/(0*)  . 

When  noise  is  present,  the  vector  forms  of  (6.26)  and  (6.27)  become 

M  +  Aa  r  (6-41) 

a/  =  a ./  +  Aa /  .  (6.42) 

Let  A L  and  A H  be  the  errors  in  the  estimation  of  the  ship’s  length  L  and  draft  H. 
It  can  be  proved  that  A L  and  A  H  will  cause  errors  AW r  and  AW/  in  the  matrices 
Wr  and  W/,  and  the  matrices  with  errors  can  be  expressed  as 


w  R  =  Wr  +  AWr  (6.43) 

W/  =  W/  +  AW/  .  (6.44) 

For  the  case  where  noise  or  errors  are  present,  (6.24)  and  (6.25)  now  become 

=  Wrc  (6.45) 

a/  =  W/d  .  (6.46) 


126 


where  the  unknowns  c  and  d  are  also  written  into  two  parts,  their  true  values  and 
errors,  i.e.,  c  =  c-f  Ac,  and  d  =  d  +  Ad.  Substitution  of  (6.41)  ~  (6.44)  into  (6.45) 
and  (6.46)  yields  two  expressions 

Ac  =  W*1  [Aafi- AWrc]  (6.47) 

Ad  =  W71  [Aa/- AW/d]  .  (6.48) 

Equations  (6.47)  and  (6.48)  show  that  the  errors  of  the  solution  of  (6.45)  and  (6.46), 
Ac  and  Ad,  depend  not  only  directly  on  the  errors  of  a/*,  a/,  W r  and  W j,  but  also 
on  the  properties  of  the  inverse  of  the  matrices  W r  and  W /.  When  the  matrices 
are  singular  or  ill-conditioned,  the  error  in  the  solution  will  be  very  large  even  for  a 
very  very  small  error  in  the  data  an,  a/,  L  and  H. 

6.4  Constrained  Linear  Inversion 

For  ill-conditioned  linear  problems,  neither  the  direct  inversion  nor  the  conven¬ 
tional  least  square  methods  work  well  [46].  A  more  powerful  method  to  treat  this 
kind  of  problem  is  the  method  of  constrained  linear  inversion  (CLI),  as  presented  in 
[46].  Before  discussing  how  to  apply  this  method  to  the  ship  hull  inversion  problem, 
the  concept  of  CLI  is  briefly  stated  below. 

For  a  general  one  dimensional  continuous  inverse  problem, 

y{0)  =  j  Ke(0,a)x(a)da  (6.49) 

Ja 

its  discrete  form  is  written  as 

y  =  Kex  (6.50) 

where  Ke  is  the  kernel  matrix  with  dimensions  mxn,y  =  [yi,  y2,  ...ym]  is  m  measure¬ 
ments  of  y(0),  and  x  =  [ii,  zj,  ...xn]  is  the  n  discrete  values  of  the  desired  function 
x(a). 
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To  solve  the  problem,  a  non- negative  scalar  measure,  ?(x),  of  the  deviations 
from  smoothness  in  x  is  introduced.  If  x  is  varied  until  q(x)  becomes  minimal, 
the  resulting  x  may  be  completely  smooth  in  the  sense  that  q(x)  will  be  zero. 
Most  measures  of  non-smoothness  are  simple  quadratic  combinations  of  X{.  For 
example,  q  =  —  x,)2  is  the  sum  of  the  squares  of  the  first  differences, 

q  =  Trial  xi+2  +  2x,+i  —  x,)2  is  the  sum  of  the  squares  of  the  second  differences, 
and  q  =  H"=i3(x«+3  —  3x<+j  +  3xt+i  —  x,)2  is  the  sum  of  the  squares  of  the  third  differ¬ 
ences.  These  summations  can  be  written  as  a  form  xTHx,  and  their  corresponding 
matrices  H  are  given  in  [46]  (pp.  124-127).  According  to  the  simulations  in  hull 
inverse  problems,  the  measures  of  smoothness  based  on  third  differences  are  very 
effective  in  reducing  the  effects  of  singularities  and  noise. 

The  constrained  linear  inversion  solution  is  obtained  by  minimizing  (Kex  — 
y)T(Kex  —  y)  with  a  constraint  of  q  —  xTHx,  and  it  can  be  written  as 

x  =  (K^K.  +  yH)~lK^y  (6.51) 

where  7  is  a  constraining  parameter.  Obviously,  7  =  0  leads  to  a  conventional  least 
square  solution.  In  a  broad  sense,  thus,  CLI  can  be  counted  as  a  least  square  method 
using  constraints.  The  most  appropriate  value  of  7  can  be  determined  by  computing 
the  residual  |Kex  —  y|.  If  the  residual  is  appreciably  larger  than  the  overall  error  in 
y,  then  7  is  too  large,  and  the  solution  has  been  over  constrained;  if  the  residual  is 
smaller  than  the  estimated  error  in  y,  the  solution  has  been  underconstrained. 

Returning  to  the  ship  hull  problem  described  in  (6.24)  and  (6.25),  we  note  that 
the  vectors  to  be  determined,  c  and  d,  are  the  coefficients  of  the  summation  of  basis 
functions  instead  of  the  hull  surface  itself.  Thus,  one  possible  method  is  to  select  the 
constraints  directly  based  on  c  and  d.  In  many  situations,  however,  these  constraints 
may  not  be  easy  to  select  for  lack  of  a  priori  information  about  the  vectors  c  and 
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d.  Usually,  few  properties  of  the  coefficients  are  known  in  advance  unlike  the  hull 
surface  function  they  represent.  For  this  reason,  the  smoothness  measure  can  be 
selected  based  on  the  hull  surface  function  C(z5  2)i  which  is  usually  smooth.  To  tell 
these  two  kinds  of  smoothness  constraints  apart,  the  constraints  based  on  coefficient 
vectors  are  called  the  coefficient  constraints  and  their  quantity  is  expressed  as  qc; 
the  constraints  based  on  hull  surfaces  are  called  the  surface  constraints  and  their 
quantity  is  expressed  as  qt. 

For  the  surface  constraints,  the  constrained  quantity  q,  can  be  constituted  from 
each  area  element  of  a  hull  surface  or  from  each  cut  line  on  the  hull  surface.  In  the 
following,  suppose  for  simplicity  that  q  is  constituted  from  the  measure  of  smoothness 
of  Nt  cut  lines  on  a  hull  surface  along  the  x-direction,  and  that  each  line  has  Nx 
points.  By  letting  Co  =  C(*«>*j)»  ^  h«o  =  h <>(£;,!/>),  the  discrete 

form  of  (6.32)  becomes 

Co  =  h£jC  +  h£,d  (6.52) 

and  the  values  on  the  j  th  line  of  the  hull  surface  can  be  expressed  as  a  vector 

zi  =  ICiy  b j  •••  OwF 

=  ¥ejC  +  90jd  (6.53) 

where  matrices  and  have  dimension  Nx  x  L\  and  Nx  x  L2,  respectively,  and 
are  defined  as 

*e;  =  [heij  he2j  ...  heN,j]T  (6.54) 

*e;  ±  [hel>  he2>  ...  he„JT  .  (6.55) 

With  these  definitions,  the  smoothness  measure  for  the  j  th  line  is 


?o(c,d)  =  zjHzj 


=  cTHOJd  +  dTHjjC  +  cTHcjC  +  drH«<,d  (6.56) 

where 

Haj  =  ¥j,H¥oi 
Hbj  = 

Hcj  =  ¥e; 

H*  =  *J,H¥oi  • 

Thus,  the  total  smoothness  measure  is  given  as 

JV* 

?«(c,d)  =  d) 

i=l 

=  cTHad  +  drH6c  +  crHcc  +  d^d  (6.57) 

where 

N, 

Ha  =  £H«y 
i=l 
N. 

Hk  =  £h6j 
i=i 
Nm 

Hc  =  gHci 

i=i 

N. 

Hi  =  . 

j=i 

In  general,  H  =  HT,  thus  HaJ  =  Hj  and  Ha  =  Hj.  Therefore,  if  H  =  HT,  then 

?,(c,d)  =  2crH.d  +  cTHec  +  dTHid  .  (6.58) 

Having  arrived  at  an  expression  for  q„  the  constrained  linear  inverse  problem  of 
the  ship  hull  can  be  described  as 

minimize :  (W*c  -  a*)T(W*c  -  a*)  +  (W*d  -  a,)T(W*d  -  a/) 


holding  constant:  g,(c,d) 
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The  solution  can  be  obtained  using  the  Lagrangian  multiplier  method  from 

Q( c,  d)  =  (WHc  -  a«)T(W*c  -  slr)  +  (W*d  -  a,)T(W*d  -  a/) 

+7?,(c,  d) .  (6.59) 

By  the  rules  of  matrix  calculus,  taking  =  0  and  =  0  yields 

W£W*+7HC  7Ha 

7Ha  WfW/  +  7H, 

For  a  ship  hull  symmetric  in  its  longitudinal  distribution  of  volume,  c  =  0  and 
a r  =  0,  and  thus  (6.60)  becomes 

(WjW/  +  7He,)d  =  Wj'a/  .  (6.61) 

In  many  real  situations,  the  surface  constraint  method  may  be  used  together  with 
the  coefficient  constraint  method  to  solve  the  linear  inverse  problem.  For  example, 
the  coefficient  constraint  method  is  used  to  solve  (6.60)  or  (6.61)  instead  of  using  the 
direct  inverse  method. 

6.5  Application  of  Bayes  Theorem  to  Inverse  Problems 

In  inverse  problems,  the  observation  data  contain  random  noise  because  they 
come  from  direct  or  indirect  measurements,  and  thus  it  is  reasonable  to  consider  the 
data  as  random  variables.  In  some  cases,  the  parameters  to  be  predicted  from  the 
observation  data  are  also  random  variables.  Hence,  we  wish  to  apply  the  statistic 
estimation  theory  to  inverse  problems  to  obtain  an  optimum  estimation  with  regard 
to  some  criterion. 

Now,  consider  a  general  problem 

y  =  f(x)  (6.62) 


c 

W£a* 

d 

Wja, 

(6.60) 


131 


where  y  is  a  vector  representing  the  observed  or  measured  data,  and  f(x)  is  a  vector 
depending  on  the  vector  x,  which  is  to  be  estimated  from  y.  It  is  assumed  that  the 
functional  relation  f  is  known  and  that  y  and  x  are  random  vectors.  If  x  denotes  an 
estimate  of  x  based  on  the  observation  vector  y,  then  the  estimation  error  is  defined 
as 

e(y)  =  x(y)  -  x  (6.63) 

In  the  Bayes  estimation  theory,  a  scale  function  of  x  and  x,  called  the  cost  function, 
is  defined,  and  in  many  cases  it  is  assumed  to  depend  only  on  the  error  of  estimate 
c.  For  example,  one  cost  function  which  is  used  frequently  represents  the  sum  of  the 
square  of  each  error  component  and  is  written  as 

<?<(«(y)j  =  «T(yMy)  (6.64) 

Another  cost  function  assigns  zero  cost  to  all  errors  less  than  and  assigns  a 
uniform  value  to  all  errors  larger  than  ±y,  that  is, 

GMy)l  = 

The  Bayes  estimation  is  based  on  the  rule  that  on  the  average  cost  is  as  small  as 
possible.  That  is,  the  Bayes  estimate  is  the  estimate  that  minimizes 

E{CtW  y)]}=/  /  C*Ny)Mx,y)dxdy  (6.66) 

where  £(•)  denotes  expectation,  and  p(x,y)  is  the  joint  probability  density  function 
of  x  and  y. 

For  the  above  two  cost  functions,  the  estimates  of  x  have  been  computed  in 
[44].  The  minimum  mean- square  error  (MMSE)  estimate  corresponding  to  the  cost 


i  Ky)l  >  f 


(6.65) 


o  Ky)l  <  f 
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function  in  (6.64)  is  given  as 


x(y)  =  £{x/y} 


(6.67) 


which  indicates  that  the  MMSE  estimate  is  the  conditional  mean  of  x  given  y.  The 
maximum  a  posteriori  (MAP)  estimate  corresponding  to  the  cost  function  in  (6.65) 
is  given  as  the  solution  of 


din  p(x/y) 
dx 


(6.68) 


or 


dlnpj y/x)  ,  dlnpjx)  . 

dx  x=*  dx  x~* 


(6.69) 


In  the  MAP  estimation,  the  estimate  is  the  values  of  x  at  which  the  a  posteriori 
density  p(x/y)  has  its  maximum. 

In  many  cases  of  interest  the  MAP  and  MMSE  estimates  are  equal  and  are 
optimum,  in  particular,  for  the  Gaussian  a  posteriori  density.  If  we  know  nothing 
about  x  other  than  the  values  of  y,  then  p(x)  is  a  constant,  and  an  estimate  x,  called 
the  maximum  likelihood  (ML)  estimate,  can  be  found  from  the  likelihood  equation: 

0.  (6.70) 


The  ML  estimate  corresponds  mathematically  to  the  limiting  case  of  a  MAP  es¬ 
timate  in  which  the  a  priori  knowledge  approaches  zero.  In  the  cases  where  the 
unknown  vector  x  is  not  a  random  vector,  the  estimate  x  is  also  obtained  from  the 
ML  estimation  in  which  the  likelihood  function  is  maximum. 

In  our  inverse  problems,  the  unknown  coefficient  vector  may  be  treat  as  a  non- 
random  vector;  thus,  the  ML  estimation  will  be  applied  to  the  ship  hull  problems.  In 
the  following,  it  is  assumed  that  the  observation  vector  y  has  a  multivariate  Gaussian 
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distribution.  The  Gaussian  conditional  probability  distribution  of  y  given  x  has  the 
form: 

PivM  =  (2Tls|i  exp{-i[y  -  f(x)JTS-‘ [y  -  f(x)J }  (6.71) 

where  n  is  the  dimension  of  y,  and  |S|  denotes  the  determinant  of  the  nxn  covariance 
matrix  S,  which  is  defined  as 

S  =  £{(y  -  £[y])(y  -  £[y]f  )  ■  (6.72) 

Substituting  (6.71)  into  (6.70)  yields 

I;(y  -  f(x)f  s-'(y  -  f(x))  =  0  (6.73) 

If  a  constraint  g(x)  is  included,  then  the  estimate  will  be  found  by  minimizing 

Q(x)  «  (y  -  f(x))TS-1(y  -  f (x))  +  7  g(x)  (6.74) 

The  above  results  will  be  used  in  the  linear  ship  hull  inverse  problem  described 
by  (6.24)  and  (6.25)  and  also  in  the  non-linear  ship  hull  inverse  problem  which  will 
be  discussed  in  the  next  subsection.  We  introduce  the  following  vectors  and  matrices 
for  the  problem  given  in  (6.24)  and  (6.25): 


Hc  H0 
H0  Hd 


With  the  above  definitions,  the  problem  in  (6.24)  and  (6.25)  is  described  as 


y  =  f  (x)  =  Wax 


(6.75) 
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and  the  constraint  in  (6.25)  is  described  as 

q(x.)  =  xtH*x  .  (6.76) 

Substituting  (6.75)  and  (6.76)  into  (6.74)  and  minimizing  the  resulting  Q(x)  yields 

(WfS~1Wa  +  7H,)x  =  WfS-'y  .  (6.77) 

Finally,  the  estimate  x  can  be  obtained  by  solving  this  equation.  Specifically,  when 
the  measurement  noise  in  the  components  of  y  are  uncorrelated  and  have  the  same 
variance  <72,  the  covariance  matrix  is  given  as  S  =  a2I,  where  I  :s  an  identity  matrix, 
and  (6.77)  becomes 

(WjW.  +  V H,)*  =  Wjy  (6.78) 

where  7'  =  <t27.  Equation  (6.78)  is  the  same  as  (6.60)  except  that  7  in  (6.60)  has 
been  replaced  with  7\  This  difference  indicates  that  the  choice  of  the  constraining 
parameter  should  include  the  consideration  of  the  noise  variance.  If  large  noise  is 
involved  in  data,  large  constraining  parameter  7‘  is  required. 

6.6  Non-Linear  inversion 

In  the  above  linear  inverse  problem,  it  has  been  assumed  that  both  the  real 
and  imaginary  parts  or,  equivalently,  both  the  phase  and  magnitude  of  the  wave 
amplitude  function  A(6)  are  well  known.  In  real  situations,  however,  the  phase 
information  of  the  wave  amplitude  function  may  not  be  available,  and  only  the 
magnitude  |>4(0)|  =  ( A7R(6)  +  A2(0)]a  is  given.  For  example,  the  phase  information 
may  be  lost  or  can  not  be  recovered  fully  when  ship  wave  spectra  are  transformed 
from  remotely  sensed  images,  or  when  wave  spectra  are  transformed  from  wave  height 
data  without  knowing  the  ship  position.  If  the  phase  information  is  not  available, 


135 


then  the  inverse  problem  becomes  more  complicated  since  it  involves  a  non-linear 
problem. 

The  non-linear  inverse  problem  can  be  derived  from  the  previously  discussed 
linear  problem.  By  squaring  both  sides  of  (6.24)  and  (6.25)  for  each  k  (k  =  1, ...,  K) 
and  then  adding  the  squared  values  for  each  k  together,  a  set  of  non-linear  equations 
is  obtained  in  the  form 


where 


A'm(9k)  =  crWdfec  +  d^W^d 
k-  1,2,3, 

A'm(9k)  £  A'Z(0k)  +  A'?(9k) 

W c*  = 

W  m  =  w/tw^ 


(6.79) 


A 

A 

To  write  the  above  equations  in  a  concise  form,  let 

A'je  t) 

a  a'm) 

— 

.  A'm(eK) 

wcfc 

wdfe 
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The  problem  described  in  (6.79)  then  becomes 

xTWmlx 


Hm  =  f(x)  = 


xTWm2x 


(6.80) 


xTWm*x 

Equation  (6.80)  is  a  general  form  for  the  non-linear  hull  surface  inversion.  If  the 
ship  hull  surface  function  can  be  written  in  the  form  of  separation  of  variables  as  in 
(5.69),  then  the  dimension  of  the  vector  to  be  solved  can  be  reduced  from  AfxJVto 
M  +  N.  This  will  result  in  a  great  reduction  of  computational  efforts,  especially  for 
large  M  and  N.  Under  the  assumption  of  separation  of  variables,  the  hull  surface 


function  ((x,  z)  can  be  expressed  in  the  form 


M  N 

C(*,2)  =  [Ea«Mx)][Ea*i&(z)] 

*=i  i=l 


(6.81) 


where  a*,-  and  atj  are  the  coefficients  to  be  determined,  and  <i>i(x)  and  flj(z)  are  the 
basis  functions.  Substitution  of  (6.81)  into  (6.1)  gives  the  equations 

(6.82) 


A'r(0)  =  t  E  Wxm{0)  ][£;«*  WZJ(0)  ] 

i=i  j=i 

Ai(0)  =  (!>*  Wxii(0) ]  [ £ WZj(0) ] 

•=i  ;=i 


(6.83) 


where  WxRi,  Wxn  and  Wzj  have  been  defined  in  (6.8)-(6.10).  For  K  values  of  0, 
there  is  a  set  of  equations 

4*(0*)  =  (xf  w***)  (x^w  zk)  (6-84) 

A'i{0k)  =  (xf  wx/fc)  (x^w  zk)  (6.85) 

k  =  1,2,---, A" 


where 


Xj  —  [q*1>  Q*2>  •••OtAf] 
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Xj"  —  [o*li  fl*2j 

wxr*  =  WxRiiPk),  Wxm{9k),  — » W*KA/(0k)]r 
w Xlk  =  [WOf/x(^/fe),  W0f/2(^fc),  •••>  WxiAf(@k)]T 
Wzk  =  [^2i(«fc),^z2(^),...,^(«fc)]T  .  (6.86) 

By  squaring  both  sides  of  (6.84)  and  (6.85)  for  each  k  (k  =  1, ...,  K)  and  then  adding 
the  squared  values  for  each  k  together,  a  vector  form  of  the  non-linear  equations  is 
obtained  in  the  form 

xf  WxiXixI’W^Xa 
xfWx2XiX^Wz2X2 

xf  W XK  XlX^  W ZK  *2 
where  W Xk  is  an  M  x  M  matrix  given  as 

w Xk  =  ^XRk  +  wXIk  wXIk 

and  Wzk  is  an  N  x  TV  matrix  given  as 

Wz*  =  wzfcw|*  . 

The  vector  equation  in  (6.87)  can  be  also  written  in  the  form 

xTW0iXXTW&ix 

xrWo2xxTWMx 

xTWatfXXTWwfX 


W  Xk  0  MN 
0  mn  ®nn 


where 
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and 

A  OMM  0  mn 

Wfcfc  = 

Omn  W Zk 

In  the  above  equations,  the  subscripts  of  0  denote  the  dimension  of  the  null  matrices. 

So  far,  we  have  established  the  non-linear  equations  for  hull  surface  inversion; 
the  next  step  we  need  to  do  is  to  solve  the  equations  to  obtain  x.  Using  the  method 
of  ML  with  a  constraint  q(x),  the  estimate  x  can  be  found  by  solving  the  following 
extremum  problem: 

minimize  :  (a^  -  f(x))rS~1(am  -  f(x)) 

holding  constant :  q(x)  .  (6.89) 

To  solve  this  extremum  problem,  several  methods  of  non-linear  optimization  can  be 
used. 

Before  we  give  examples  of  ship  hull  estimations,  it  is  worth  to  have  some  discus¬ 
sions  about  the  issue  of  uniqueness.  Like  many  other  inverse  problems,  the  unique¬ 
ness  problem  arises  in  the  inversion  of  a  ship’s  hull  shape,  that  is,  whether  or  not  a 
hull  form  can  be  uniquely  determined  form  a  free  wave  spectrum.  We  may  see  this 
problem  from  two  aspects,  that  is,  from  physical  and  mathematical  aspects. 

In  physics,  the  question  is  whether  different  hull  shapes  have  their  respective 
wave  patterns.  According  to  Newman’s  study  on  this  problem  [9],  there  is  one  and 
only  one  equivalent  source  distribution  for  a  given  wave  system;  however,  different 
physical  source  distributions  can  be  related  to  the  same  equivalent  source  distribu¬ 
tion,  different  vessels  can  be  responsible  for  the  same  wave  system,  and  hence,  in 
general,  non-uniqueness  exists  in  the  hull  inversion.  In  mathematical  processing,  the 
uniqueness  problem  may  also  arises  even  if  an  original  inverse  problem  is  unique. 
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For  example,  the  non-linear  inversion  of  hull  shapes  described  above  may  have  mul¬ 
tiple  solutions  due  to  the  non-linear  operations;  thus  it  is  more  complicated  than  the 
above  linear  inverse  problem  of  hull  shapes. 

Because  of  the  above  reasons,  we  must  apply  proper  constraints  on  hull  form 
to  improve  problem  conditioning.  Such  constraints  may  facilitate  a  unique  relation 
between  a  free  wave  spectrum  and  a  hull  form  [8].  Although  we  have  no  idea  about 
a  hull’s  exact  shape  before  we  do  the  inversion  for  it,  we  may  have  some  useful 
knowledge  about  man-made  ship  hulls  and  use  it  as  constraints.  For  example,  hull 
surface  smoothness  constraints,  closure  type  constraints  and  volume  constraints  axe 
useful  in  hull  inversion.  In  the  following  simulations,  the  hull  surface  smoothness 
constraint,  especially  the  third  difference  smoothness  constraint,  is  used. 

6.7  Examples  of  Ship  Hull  Estimations 


The  previous  sections  have  discussed  the  problems  of  a  ship  hull’s  linear  and 
non-linear  inversion  and  the  methods  to  solve  them.  This  section  evaluates  the  hull 
model  given  in  (6.32),  and  gives  the  simulation  results  of  mathematically  well  defined 
hulls  amd  the  Quapaw  hull. 

To  evaluate  the  error  performance  of  the  estimated  hull  surface  £(x,z),  we  con¬ 
sider  the  absolute  error  ea  =  |Cy  —  C<;|»  the  relative  error  e,  =  |C.j  —  C«jl/IC«j|i  and  the 
relative  overall  r.m.s.  residual  error  which  is  defined  as 


A 

^ro  — 


f  J[Z(x,z)-C{x,z)]2dxdz 
\  f  f  C2{x,z)dxdz 


(6.90) 


or,  in  the  discrete  form,  as 


^rn  — 


I  Cij)2 

N  XijCb 


(6.91) 


There  is  some  difficulty  in  using  the  relative  error  er  to  evaluate  the  performance  on 


the  edges  of  the  hull  surface  because  the  true  values  on  the  edges  approach  zero.  The 
relative  overall  r.m.s.  residual  error  can  be  envisaged  as  the  ratio  of  the  volume  of 
the  error  to  the  volume  of  the  hull. 

In  addition  to  the  above  criterion  on  hull  surface  error  performance,  we  may  also 
consider  the  residual  error  due  to  the  inversion  operation: 

€re,  =  (a R  -  W„  c)T(a* -W„c)  +  (a/  -  W,  d)r(aj  -  W,  d)  .  (6.92) 

Note  that  a  large  error  of  cre,  results  in  a  large  error  of  ero.  However,  even  a  very 
small  error  of  ere<  never  guarantees  a  small  error  of  e™,  that  is,  e™  may  be  very  large 
even  for  a  very  small  erM. 

6.7.1  Evaluation  of  Hull  Surface  Models 

Before  giving  the  simulation  results,  we  first  evaluate  the  performance  of  the  hull 
surface  model.  The  basis  functions  selected  from  (6.17)  and  (6.20)  consist  of  only 
polynomials;  thus  the  hulls  defined  with  polynomials  can  be  exactly  expressed  when 
the  numbers  of  basis  functions  are  sufficiently  large.  So,  we  do  not  consider  these 
types  of  hulls  here,  but  the  Wigley-Cosine  hull  and  the  Quapaw  hull. 

For  the  evaluation  of  a  hull  surface  model,  the  model  vector  is  calculated  from 
(6.32)  using  the  least  square  method  given  the  values  of  £(x,  z),  then  the  resulting 
model  vector  xm  =  [cm  dm]r  is  used  to  evaluate  the  hull  surface  values,  (m(x,  z)  = 
(x,  z)cm  +  hf  (x,  z)dm. 

The  Wigley-Cosine  hull  has  length  L  =  100  meters,  draft  H  —  10  meters  and 
width  B  =  10  meters.  Table  6.1  lists  the  values  of  the  relative  overall  r.m.s.  residual 
error  ero  for  different  numbers  of  basis  functions.  The  relative  overall  r.m.s.  residual 
error  ero  versus  the  number  of  basis  functions  N  is  shown  in  Figure  6.4,  assuming 
that  the  numbers  of  basis  functions  in  the  x-  and  z-directions  are  equal,  i.e.,  M  =  N. 
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€ro(%) 

N=1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

M  =  1 

48.82 

21.95 

18.14 

18.14 

18.14 

18.14 

18.14 

18.14 

18.14 

18.14 

M  =  2 

47.18 

16.88 

11.37 

11.37 

11.37 

11.37 

11.37 

11.37 

11.37 

11.37 

M  =  3 

46.23 

13.16 

3.96 

3.96 

3.96 

3.96 

3.96 

3.96 

3.96 

3.96 

S 

II 

*»■ 

46.12 

12.66 

1.61 

1.61 

1.61 

1.61 

1.61 

1.61 

1.61 

1.61 

K 

II 

Oi 

46.11 

12.63 

1.33 

1.33 

1.33 

1.33 

1.33 

1.33 

1.33 

1.33 

M  =  6 

46.10 

12.57 

0.51 

0.51 

0.51 

0.51 

0.51 

0.51 

0.51 

0.51 

M  =  7 

46.10 

12.57 

0.51 

0.51 

0.51 

0.51 

0.51 

0.51 

0.51 

0.51 

M  =  8 

46.09 

12.56 

0.26 

0.26 

0.26 

0.26 

0.26 

0.26 

0.26 

0.26 

M  =  9 

46.09 

12.56 

0.26 

0.26 

0.26 

0.26 

0.26 

0.26 

0.26 

0.26 

M  =10 

46.09 

12.56 

0.15 

0.15 

0.15 

0.15 

0.15 

0.15 

0.15 

0.15 

Table  6.1:  The  Wigley-Cosine  hull’s  relative  overall  r.m.s.  residual  error  ero  in  the 
hull  approximation  using  M  x  N  basis  functions. 
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tro(%) 

N=1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

M  =  1 

45.62 

26.16 

21.77 

20.67 

20.50 

20.50 

20.40 

20.28 

20.22 

20.21 

M  =  2 

44.46 

22.58 

17.16 

15.34 

14.62 

14.41 

14.32 

14.31 

14.24 

14.23 

M  =  3 

44.33 

20.68 

11.87 

8.99 

7.77 

7.26 

7.02 

6.82 

6.69 

6.63 

II 

S 

44.32 

20.57 

11.04 

7.89 

6.42 

5.75 

5.44 

5.33 

5.23 

5.13 

M  =  5 

44.21 

20.35 

10.58 

6.85 

5.02 

4.07 

3.61 

3.37 

3.28 

3.17 

M  =  6 

44.15 

20.20 

10.32 

6.27 

4.20 

2.99 

2.33 

1.99 

1.85 

1.69 

II 

S 

44.14 

20.20 

10.29 

6.25 

4.06 

2.78 

2.04 

1.67 

1.45 

1.40 

s 

II 

00 

44.13 

20.19 

10.28 

6.24 

3.99 

2.67 

1.89 

1.50 

1.25 

1.19 

M  =  9 

44.13 

20.19 

10.27 

6.22 

3.97 

2.61 

1.81 

1.38 

1.16 

1.05 

M  =10 

44.11 

20.16 

10.24 

'  6.17 

3.91 

2.51 

1.66 

1.17 

0.91 

0.78 

Table  6.2:  The  Quapaw’s  relative  overall  r.m.s.  residual  error  Cr0  in  the  hull  approx¬ 
imation  using  M  x  N  basis  functions. 
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Figure  6.4:  The  relative  overall  r.m.s.  residual  error  of  the  Wigley-Cosine  hull 
and  the  Quapaw  model  in  the  hull  approximation  using  N  x  N  basis 
functions. 

The  Quapaw  hull  has  length  L  =  4.953  meters,  draft  H  =  0.362  meters  and  width 
B  =  0.978  meters.  Table  6.2  lists  the  values  of  the  relative  overall  r.m.s.  residual 
error  er0  for  different  numbers  of  basis  functions.  The  curve  of  eT0  versus  N  with  the 
condition  of  M  —  N  is  also  plotted  in  Figure  6.4. 

From  the  above  analysis,  the  hulls  can  be  approximately  very  well  for  the  suffi¬ 
ciently  large  numbers  of  the  basis  functions.  For  example,  e™  can  be  less  than  3.0% 
for  M  =  N  =  6,  and  less  than  1.0%  for  M  =  N  =  10.  In  actual  situations,  the 
noise  in  the  wave  amplitude  function  may  result  in  larger  error  in  the  recovered  hull 
surface  than  this  caused  by  the  model  approximation. 
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6.7.2  Estimation  of  Hull  Shape  from  the  Wave  Amplitude  Function 

It  is  assumed  in  this  subsection  that  both  the  real  and  imaginary  parts  of  the 
wave  amplitude  function  are  fully  available.  Thus,  the  model  vectors  c  and  d  are 
estimated  from  the  linear  inversion. 

In  the  following,  we  consider  first  a  hull  defined  with  polynomials: 

C(x,z)  =  y  (1  -x  -  xJ  +  x3)(l  -  z2)  x  €  [-1,1],  z  €  [-1,0]  (6.93) 

which  is  more  typical  than  Wigley’s  parabolic  hull.  This  hull  is  not  symmetrical  in 
the  longitudinal  direction,  and  thus  its  wave  amplitude  function  contains  both  real 
and  imaginary  parts. 

In  the  simulation,  the  values  of  the  wave  amplitude  function  A{0)  are  calculated 
from  (5.8)  with  wave  angles  0  in  [5°,  80°].  The  basis  functions  are  selected  from(6.17) 
and  (6.20).  With  the  basis  functions  and  the  given  ship  length  and  draft,  the  kernel 
matrices  W«  and  W /  can  be  computed.  Once  the  hull  model  vector  x  =  [c  d]T  is 
determined  by  solving  the  linear  inversion  problem,  the  hull  surface  values  can  be 
evaluated  from  (6.32),  i.e.,  £(x,  z)  —  h^(x,  z)c  +  h^(x,  z)d.  In  the  calculations  below, 
all  the  constraints  are  based  on  the  measures  of  the  sum  of  the  squares  of  the  third 
differences. 

With  the  above  procedures,  simulation  results  can  be  obtained.  As  an  example, 
it  is  assumed  that  the  ship  speed  is  U  =  10  meters/second,  the  parameters  in  (6.93) 
are  L  =  100  meters,  H  =  10  meters,  and  B  =  10  meters.  The  dotted  line  shown  in 
Figure  6.6  is  the  wave  amplitude  function  calculated  from  (5.8). 

When  the  numbers  of  basis  functions  are  taken  as  M  =  4  and  N  =  3,  the 
constrained  ML  method  gives  an  exact  solution  for  the  containing  parameters  7r  =  0 
and  7,  =  10-9.  As  the  numbers  of  basis  functions  increase,  the  error  appears  in  the 
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solution,  because  the  kernel  matrices  become  singular.  Figure  6.5  shows  the  curves 
of  ero  and  ere,  versus  7,  when  M  =  6  and  N  =  6.  It  can  be  seen  from  the  figure  that 
there  is  a  minimum  value  1.0557%  of  e™  based  on  the  calculation  of  11  x  11  points 
on  the  hull  surface,  associated  with  a  residual  error  =  7.4  x  10-8.  The  estimated 
hull  surface  is  plotted  in  Figure  6.7  together  with  the  error  surface  contour  plots  of 
ea  and  e,. 

To  observe  the  effects  of  the  error  in  the  hull  length  and  draft  and  the  noise  in 
the  wave  amplitude  function  on  the  estimation  performance,  a  couple  of  assumptions 
are  made:  first  that  there  is  10%  error  in  the  ship  length  and  draft,  that  is,  L  —  110 
and  H  =  11  meters  are  used  in  the  computation  of  the  kernel  matrices;  second,  that 
there  is  10%  noise  in  the  real  and  imaginary  parts  of  the  wave  amplitude  function 
based  on  the  energy,  that  is,  a  10  dB  signal  to  noise  ratio.  When  the  noise  is  assumed 
to  have  a  Gaussian  distribution  with  zero  mean,  the  standard  deviation  is  taken  as 


<rnR  —  yO.l (o\R  +  E[Ar\2)  for  the  real  part,  and  <rnj  =  y0.1(<r^7  +  E[Ai]2)  for  the 
imaginary  part  of  the  wave  amplitude  function.  The  wave  amplitude  function  with 
noise  is  shown  with  the  solid  lines  in  Figure  6.6.  Figure  6.8  shows  the  estimated  hull 
surface  together  with  the  error  surface  contour  plots  of  ea  and  c,,  given  M  =  N  =  6. 
In  the  calculation,  both  the  coefficient  constraint  and  surface  constraint  are  used, 
with  7C  =  10-1  and  7,  =  380.  The  overall  error  performances  are  given  as  eT0  =  27.2% 
and  cre»  =  0.482. 


With  the  same  parameters  and  conditions,  we  again  do  the  simulation  for  the 
Wigley-Cosine  Hull.  The  wave  amplitude  function  and  the  simulation  results  are 
displayed  in  Figure  6.9  and  Figure  6.10.  The  overall  error  performances  are  given  as 
Cro  =  32.5%  and  cre,  =  0.55,  with  7C  =  10~2  and  7,  =  160.  Figure  6.11  shows  the 
curves  of  tro  and  tre,  verses  7,,  given  7C  =  10-1  and  10“ 2,  for  the  modified  Wigley 
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Figure  6.5:  Curves  of  ero  and  ere«  versus  7,  for  the  modified  Wigley  hull  with  M  = 
N  =  6.  No  error  and  noise  are  present  in  the  input  data. 

hull  described  in  Eq.(6.93)  and  the  Wigley-Cosine  Hull.  It  is  found  from  the  figure 
that  the  error  of  ere4  has  only  a  small  change  for  different  7,  with  a  given  7C.  However, 
has  obvious  variations  for  different  je. 
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Figure  6.6:  (a)  the  real  part,  (b)  imaginary  part,  and  (c)  the  magnitude  of  the  wave 
amplitude  function  of  the  modified  Wigley  hull,  given  the  parameters 
L  =  100  m,  H  =  10  m,  B  =  10  m,  and  U  —  10  m/s.  The  dotted  lines 
represent  the  wave  amplitude  function  with  no  noise,  and  solid  lines  the 
one  with  10%  noise. 


z(m) 


W  True  Hun  Surface 


VeneCs 


0>)  Estimated  HbJJ  Surface 


Vo**C\ 


4fe2jgje  error  (W 


7 

~50  -40  -30  .20 


.o.o\- 

.0.01 


*20  o 

X  (m) 


Olj 
-2 1 
-4  ff 
-6 1 
-8 1 
-10§ 


»  — ^ 

Jglative  error 


10  20  30  40  50 


-40  -30  -20  -'vpr-r1r^0 

figure  6. 7:  The  modified  w  M 

and  7,  =  10~«  ’  ^  ~  10  m>  ^  =  10  m/s  j|/L  £  P^ameters 

’  d,  iv  =  6;  7c  - 


149 


(*)  Tree  Hull  Surf** 


<b)E*tfa,««d  Hon  Surface 


r^^h^error(m) 


W  x(m)  10  20  30  40  50 

'/mr^W?lativews(%) 


Figure  6  8- Th  ”(m)  W  20  30  40  50 

8 The  modified  Widev  h„ii 

“<l  7.  =  3800.  m’  B  =  10  m,  V  =  l0  m/s  e"  ,ke  Parameters 

’  /V  =  6>  7c  =  lO"1 


150 


(a) 


0  10  20  30  40  50  60  70  80  90 

B  (deg.) 


0  10  20  30  40  50  60  70  80  90 

B  (deg.) 


Figure  6.9:  (a)  The  real  part,  (b)  imaginary  part,  and  (c)  the  magnitude  of  the 
wave  amplitude  function  of  the  Wigley-Cosine  hull,  given  the  parameters 
L  =  100  m,  H  =  10  m,  B  =  10  m,  and  U  —  10  m/s.  The  dotted  lines 
represent  the  wave  amplitude  function  with  no  noise,  and  solid  lines  the 
one  with  10%  noise. 
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Figure  6.11:  Curves  of  ero  and  cre,  verses  7,  given  7C  =  10-3  and  10~2,  and  the 
parameters  L  =  100  m,  H  —  10  m,  B  =  10  m,  U  =  10  m/s,  M  =  N  =  6, 
a  10%  noise  or  error  in  the  input  data.  In  the  subscripts  of  e,  “mw” 
represents  the  modified  Wigley  hull,  “wc”  the  Wigley-Cosine  hull,  “a” 
denotes  7C  =  10-1,  “b"  7C  =  10-2. 
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6.7.3  Estimation  of  Hull  Shape  from  the  Magnitude  of  the  Wave  Ampli¬ 
tude  Function 

It  is  assumed  in  this  subsection  that  only  the  magnitude  of  the  wave  amplitude 
function  is  available.  Thus,  the  model  vectors  c  and  d  will  be  obtained  by  solving  the 
extremum  problem  described  in  (6.80)  and  (6.89).  In  the  following,  the  conjugate 
gradient  method  [42]  is  used  to  minimize  Q(x)  =  (Bm— f(x))TS-1(am— f(x))+79.(x), 
where  g,(x)  is  a  surface  smoothness  constraint.  This  multiple-dimensional  extremum 
problem  may  result  in  multiple  solutions,  thus  the  constraint  must  be  considered. 
Additionally,  an  initial  vector  of  x  needs  to  be  imposed  for  the  conjugate  gradient 
method.  To  improve  problem  conditioning,  the  initial  values  of  the  coefficient  vector 
x  are  obtained  using  the  method  as  shown  in  Subsection  6.7.1  from  a  known  hull 
model,  given  the  same  hull  length  and  draft  and  the  same  numbers  of  basis  functions 
which  will  be  used  in  the  inversion.  In  the  examples  below,  the  variables  of  the  hull 
surface  functions  are  assumed  not  be  separable,  the  numbers  of  basis  functions  are 
taken  as  M  =  6  and  N  =  6,  and  all  the  constraints  are  based  on  the  measures  of  the 
sum  of  the  squares  of  the  third  differences. 

The  first  example  is  the  Wigley- Cosine  hull  with  the  same  hull  parameters  as 
before  but  with  no  noise  in  the  input  data.  The  magnitude  of  the  wave  amplitude 
function  is  shown  as  a  dotted  line  in  Figure  6.9(c).  The  initial  values  of  the  coefficient 
vector  is  obtained  from  the  modified  Wigley  hull  given  in  (6.93),  and  no  constraint  is 
considered.  The  recovered  hull  is  given  in  Figure  6.12  together  with  its  error  surface 
contour  plots  of  c0  and  er.  The  overall  error  performances  are  given  as  er0  =  8.92% 
and  cre,  =  0.0568. 


The  next  example  is  the  same  Wigley-Cosine  hull,  but  10%  noise  is  considered. 
The  magnitude  of  the  wave  amplitude  function  with  noise  is  shown  as  a  solid  line  in 
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33.47 

10.0 

10.0 
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100.0 

0.8977 

28.85 

Table  6.3:  The  error  performance  of  the  Wigley-Cosine  Hull  estimated  from  the  mag¬ 
nitude  of  the  wave  amplitude  function  with  M  =  N  =  6. 

Figure  6.10(c).  The  hull  and  draft  also  have  a  10%  of  error.  The  initial  values  of  the 
coefficient  vector  is  obtained  from  the  modified  Wigley  hull  model  given  in  (6.93), 
and  the  surface  constraint  is  considered  with  7,  =  100.  The  recovered  hull  is  shown 
in  Figure  6.13  together  with  its  error  surface  contour  plots  of  ta  and  er.  The  overall 
error  performances  are  given  as  e,0  =  28.85%  and  tre,  =  0.9037.  Table  6.3  shows 
the  simulation  results  under  different  conditions,  in  which  e/,,  tjj,  Ca„  represent  the 
relative  errors  in  L,  H  and  Am(0),  respectively. 

Finally,  we  consider  the  Quapaw  hull,  a  more  practical  example  in  which  the  wave 
amplitude  function  is  calculated  from  the  spectrum  of  the  wave  elevation  measured 
from  the  tow  tank  as  seen  in  the  previous  two  chapters.  A  total  of  40  data  points  of 
the  magnitude  of  the  wave  amplitude  functions  is  used  in  the  hull  surface  estimation, 
corresponding  to  wave  angles  in  [12°,  70°].  Figure  6.14  shows  the  magnitude  of  the 
wave  amplitude  function  verse  wave  angles  9  for  the  data  RUN5-A  and  RUN5-B. 
The  hull  surface  estimated  from  RUN5-A  and  RUN5-B  are  shown  in  Figure  6.15  and 
Figure  6.16  together  with  the  error  surface  contour  plots  of  eQ  and  er  for  M  =  N  =  6. 
The  overall  error  performances  are  given  as  tTO  =  34.607%  and  erM  =  1.335  for  RUN5- 


1 


I 

I 

I 

I 


A,  and  tro  —  34.662%  and  tre,  =  1.250  for  RUN5-B.  Figure  6.17  shows  the  curves  of 
tro  and  erej  verses  7,. 

In  the  above,  we  have  demonstrated  the  examples  of  ship  hull  estimation  from 
the  wave  amplitude  functions  and  their  magnitude.  From  these  examples  and  the 
simulations  with  different  parameters  and  conditions,  some  comments  can  be  made 
as  follows. 

(1)  When  input  data  are  perfect  and  the  hull  to  be  recovered  can  be  exactly  ex¬ 
pressed  using  the  basis  functions,  the  exact  solution  of  this  hull’s  surface  ((x,  2)  can 
be  obtained  in  the  linear  inversion  with  reasonable  numbers  of  basis  functions,  for 
example,  Af  =  4  and  N  =  3  for  the  hull  given  in  (6.93).  However,  there  will  be  errors 
in  the  solution  when  the  input  data  are  not  accurate  and  have  noise. 

(2)  In  linear  inverse  problems,  the  estimation  performance  is  usually  much  more 
sensitive  to  the  noise  in  the  wave  amplitude  function  than  the  errors  on  the  hull’s 
length  and  draft.  Larger  constraining  parameters  axe  needed  to  achieve  a  better 
result  for  larger  noise.  In  addition  to  a  surface  constraint,  a  coefficient  constraint  is 
helpful  to  achieve  a  stable  solution,  especially  in  a  severe  noise  environment. 

(3)  The  smoothness  measure  based  on  the  sum  of  the  squares  of  the  third  differences 
are  more  effective  than  some  other  measures  such  as  those  based  on  the  sum  of  the 
squares  of  the  first  or  second  differences. 

(4)  The  error  in  a  hull  length  L  and  draft  H  mainly  affects  the  accuracy  of  the  side 
edge  and  bottom  edge  of  the  hull  in  solutions.  In  general,  the  relative  error  in  the 
solutions  is  larger  on  the  hull’s  bottom  boundary  than  those  on  the  upper  boundary 
and  central  part. 
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Figure  6.14:  The  magnitude  of  the  wave  amplitude  function  obtained  from  the  data 
RUN5-A  and  RUN5-B  for  the  Quapaw  model  with  L  =  4.953  m,  H  = 
0.362  m,  B  =  0.978  m,  and  U  =  2.229  m/s. 
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Figure  6.17:  Error  performance  of  the  Quapaw  hull  estimated  from  the  data  RUN5- 
A  and  RUN5-B,  given  M  =  N  =  6.  The  subscript  “a”  represents 
RUN5-A,  and  ub”  RUN5-B. 


CHAPTER  VII 


CONCLUSIONS  AND  RECOMMENDATIONS 


The  previous  chapters  have  presented  the  study  on  the  estimation  of  a  moving 
ship’s  speed,  direction,  length  and  hull  shape  from  its  wave  spectra.  Deep  water, 
steady  waves  and  the  linearized  free  surface  condition  are  assumed  throughout  the 
study. 

The  estimation  of  ship  speed  and  direction  has  relied  on  the  distinct  features  of 
the  Fourier  spectrum  of  the  ship  generated  waves.  The  introduction  of  the  concept 
of  the  complex  wave  elevation  has  been  useful  in  simplifying  the  derivation  of  the 
wave  amplitude  function  from  one  and  two  dimensional  wave  spectra.  In  general,  it 
is  suggested,  for  high  accuracy  and  the  ability  to  separate  background  noise,  that  a 
ship’s  speed  and  direction  be  estimated  from  two  dimensional  wave  spectra.  However, 
for  high  accuracy  and  easy  computation,  we  suggest  that  the  wave  amplitude  function 
be  recovered  from  one  dimensional  wave  spectra. 

The  extraction  of  a  ship’s  hull  geometry  has  been  based  on  the  relation  of  the 
wave  amplitude  function  and  the  ship  geometry  under  the  assumption  of  the  thin- 
ship  theory.  The  theoretical  model  developed  for  the  wave  amplitude  function  has 
explicitly  revealed  the  periodic  character  of  the  wave  amplitude  function  and  its  re¬ 
lation  with  the  ship  length.  From  this  model,  it  is  also  found  that  the  ship’s  bow 
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and  stern  play  a  dominant  role  in  the  wave  amplitude  function,  and  that  disconti¬ 
nuities  in  the  hull  surface  function  and  its  derivatives  results  in  the  increase  of  other 
frequency  components.  The  real  and  imaginary  parts  and  the  magnitude  of  the  wave 
amplitude  function  generally  appear  to  be  signals  which  are  both  magnitude  and 
angle  modulated.  The  bow  and  stern’s  concaveness  or  convexity  have  effects  on  the 
frequency  variation,  or  the  distribution  of  the  zero-crossing  points  of  the  wave  ampli¬ 
tude  function,  and  may  result  in  over-  or  underestimation  of  ship  length.  The  three 
methods,  the  spectrum  method,  zero-crossing  method  and  frequency  demodulation 
method,  have  shown  their  effectiveness  in  the  estimation  of  a  ship’s  length. 

The  extraction  of  a  ship’s  geometry  information  is  essentially  a  linear  or  nonlinear 
inverse  problem.  The  ill-conditioning  of  the  kernel  matrices  is  a  critical  problem  in 
the  linear  inversion.  It  has  been  shown  that  the  constrained  maximum  likelihood 
method  for  hull  surfaces  or/and  coefficients  is  useful  in  reducing  the  effects  from  this 
ill-conditioning  and  the  noise  present  in  input  data.  When  the  noise  components, 
present  in  the  data  of  the  wave  amplitude  function,  have  independent  identical  Gaus¬ 
sian  distributions,  the  constrained  maximum  likelihood  method  and  the  constrained 
linear  inversion  method  are  equivalent  in  the  case  of  linear  inversion.  The  uniqueness 
of  the  solution  is  another  critical  problem  in  the  nonlinear  inversion.  In  addition  to 
holding  constraints,  the  proper  choice  of  the  initial  coefficient  vector  has  been  shown 
to  be  helpful  in  finding  the  desired  solution.  In  the  examples  given  in  Chapter  6, 
the  initial  vectors  were  chosen  from  the  known  hull  models  given  the  parameters  of 
length  and  draft.  According  to  simulation  results,  high  accuracy  can  be  achieved 
when  the  hulls  are  recovered  from  the  perfect  input  data.  However,  the  accuracy 
decreases,  especially  at  the  edges  of  hulls,  when  noise  is  present  in  input  data. 

The  results  from  the  study  show  a  promising  possibility  of  detecting  a  moving 
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ship’s  characteristics  through  measuring  the  ship  generated  wave  pattern.  For  a  prac¬ 
tical  application  of  this  technique,  however,  many  problems  in  theory  and  practice 
still  remain  to  be  investigated  and  solved.  The  following  are  some  considerations  and 
recommendations  for  further  study. 

(1)  The  estimation  of  a  ship’s  length  and  hull  shape  given  in  the  study  is  based  on  the 
thin-ship  theory.  The  theory  itself  has  its  limitations  and  weakness.  The  example  of 
the  Quapaw  hull  has  shown  that  the  methods  developed  under  the  thin-ship  theory 
tire  also  suitable  for  fat  ships  like  the  Quapaw,  however,  further  investigation  of 
other  more  general  types  of  hulls  will  be  helpful  in  understanding  the  performance 
and  limitation  of  the  methods  which  have  developed. 

(2)  In  the  hull  inversion,  it  was  assumed  that  the  hull  draft  was  known.  However,  it 
is  still  an  open  problem  -  how  to  estimate  the  hull  draft  from  ship  wave  spectra. 

(3)  The  nonlinear  inversion  described  in  (6.88)  for  a  variable-separable  hull  function, 
having  a  relatively  small  size  of  coefficient  vector,  has  not  been  tested  and  needs  to 
be  further  evaluated. 

(4)  According  to  the  results  of  the  hull  inversion,  constraints  are  important  and 
necessary  for  obtaining  a  stable  solution  with  a  reasonable  accuracy  in  either  linear  or 
nonlinear  problems.  Thus,  we  need  to  search  further  for  an  optimal  set  of  constraints, 
if  they  exist,  to  achieve  better  performance  in  the  hull  estimation. 

(5)  In  the  above  study,  constant  ship  velocity  is  assumed,  thus  only  a  steady  ship 
wave  field  was  considered.  However,  it  may  not  be  always  true  in  practical  situations. 
The  effects  of  a  ship’s  velocity  variation  and  other  unsteady  effects  on  the  estimations 
have  not  been  studied.  To  estimate  a  ship’s  characteristics  from  an  unsteady  ship 
wave  field,  new  techniques  may  need  to  be  developed. 

(6)  In  addition  to  the  techniques  developed  in  this  study,  another  critical  problem 
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existing  in  the  application  of  remote  sensing  of  a  moving  ship  and  extracting  its 
characteristics  is  how  to  transform  successfully,  a  remotely  sensed  image  or  intensity 
spectrum  to  a  quantitative  wave  elevation  spectrum.  Monaldo  and  Lyzenga’s  efforts 
in  the  research  on  the  estimation  of  the  wave  slope-  and  height- variance  spectra  from 
SAR  imagery  [5]  may  be  helpful  in  solving  this  problem,  but  further  study  on  it  is 
still  needed. 
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APPENDIX 

WAVE  AMPLITUDE  FUNCTION  WITH  ONE 
DISCONTINUITY 


This  appendix  gives  the  integration  results  with  the  integral 

A(Kk)  =  j1J(x)e^dx 

and  the  integrand 


(1) 


F{x)  = 


<*jXJ  for  -1  <  x  <  xb 


(2) 


Ejto  bjX3  for  Xfc  <  x  <  1 
Substituting  (2)  into  (1)  yields  the  results 


MKX)  =  <9R,(tfr)cOs(l/  -<f>Rl  (v))  +  QRl(Kx, Xb)cOs(xbV  -  <f>R2{v,Xb))  (3) 
Ai(Kx)  =  Qi1{Kt)cos(i/-<t>Ii(v))  +  QIi(Kx,xb)cos(xbi'-<t>ii(v,xi> ))  (4) 


with 


QbAKt) 


Qn3(K x) 
Qh{Kx) 


{(I^Pft(ftl) -a.Pft(ft-l))]2  + 

«=o 

(£(&«9ft(ft  *)  +  <*.■  9Ki(ft-l))]2}’  (5) 

i'=0 

{  [£(&<  -  °i)  Pft  (ft  xb))  ] 2  +  [  O&i  “  a‘)  9ft  (ft  *&))  ]*}*  (6) 

•=0  i=0 

u£(feiP/.(fti)  -oiP/,(ft-i))]2  + 

i=0 
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Qh(K*) 

4>iM 
4>h{v,  *») 


1 1 Z(k  i)  +  <*i  qi,(v,  -i))  ]2}’ 

t=0 

{  [jt(bi  -  «i)P/i(vt  *»))]*  +  (£(fe<  -  a.)?/;^^))]2}* 

»=0  i=0 

t  -it  E.^o  [bigRifa  1)  +  g^«.-(y»  -1)  1 
E~o  [kPRifa  1)  ~  aiPRi(v> -1) 

tan"1  f  , 

E£o  (ft«  -  ~ai)PRi{^ X b) 

tan"1  r  ££° 11  +  a'gA(y»  ~1)  t 

lESoI^Kl)-a.P/,K-l)J 

tan-1  f  ^,=0  g>)  i 

E£o(&< - a.Op/^^Xfc) 


(7) 

(8) 
(9) 

(10) 

(11) 

(12) 


where  pni(a,x),  <^(0,1),  p/^a,!)  and  qi^a^x)  are  as  defined  in  (5.16)~  (5.19). 
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